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TWO-DIMENSIONAL PROBLEMS OF THE 
DECAY OF MAGNETIC FIELDS IN 
MAGNETOHYDRODYNAMICS 


By T. G. COWLING and A. HARE (The University, Leeds) 
Received 24 January 1957] 


SUMMARY 


The normal modes of decay of a magnetic field in the presence of steady motions 
in a fluid conductor are studied in the two-dimensional case. Slow motions are 
considered by perturbation methods, which show that to a first approximation the 
motion always increases the rate of decay in the lowest mode. For fast motion, 
there exists a limited class of modes little affected by the motion, in which the lines 
of force and streamlines nearly coincide. In the general class of modes the motion 
profoundly affects the decay ; these are studied in two special cases, those of stream- 
ing in a fluid slab between parallel walls and of non-uniform rotation in a circular 
cylinder. The motion is found to transport and distort the field, and to increase 
the rate of decay in the lower modes, roughly proportional to the % power of the 
velocity. For the slab problem, velocities of intermediate magnitude are studied 
by numerical methods. The extension of the results to more general classes of tavo- 
dimensional motion is briefly considered. 


1. Statement of the problem 
RECENTLY Bullard and Gellman(1) have established, by numerical methods, 
that a steady magnetic field can be maintained by dynamo interaction with 
asteady motion of suitable size and type in a conducting fluid. The details 
of the mechanism are, however, somewhat obscure. For this reason it 
appears desirable to study the normal modes of decay of magnetic fields in 
the presence of general classes of steady motions in conducting fluids. By 
this is meant modes in which the field H after time ¢ is connected with the 
initial field H, by an equation 

Hn = He, (1) 
where c may be either real or complex. Dynamo maintenance of a steady 
field corresponds to ¢ = 0; more generally, whenever the real part of o 
passes from positive to negative values as the velocity is varied, the field 
passes from a decreasing mode to an amplifying mode. 

The simplest problem to discuss is a two-dimensional one, in which 
currents flow parallel to Oz in a cylinder of uniform liquid with generators 
parallel to Oz; the velocity v and the magnetic field H are perpendicular 
to Oz. In this case dynamo maintenance is known to be impossible (2), so 
that the results have only an indirect bearing on the dynamo problem. 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 4 (1957)] 
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Nevertheless, because of that indirect bearing, and because illustrative 
results are easier to obtain than in the three-dimensional case, the two- 
dimensional problem will be considered here. Chandrasekhar has recently 
considered the three-dimensional problem with axial symmetry (3); the 
topological similarity of this with the two-dimensional problem should lead 
to similar results in the two cases. However, Chandrasekhar found that a 
sufficiently fast motion might slow down the decay almost indefinitely, 
The present work gives no support to such a view; it indicates that a rapid 
motion always leads to rapid decay. The discrepancy between the results 
appears to be due to the inadequacy of certain approximations made by 
Chandrasekhar. It can in fact be shown by general methods that no motion 
whatsoever can, either in the two-dimensional or the axisymmetric case, 
increase the time of decay beyond a definite upper limit. 

In the two-dimensional case, to avoid having an infinite total magnetic 
energy per unit length in the z-direction (which would make any decay 
problem meaningless) it is necessary to suppose that zero total electric 
current flows across any section z = constant of the conducting cylinder. 
Appropriate methods of approximation to the fields are available when the 
motion is either slow or fast. These do not apply to intermediate values 
of the velocity; to link up the results for slow and fast motions, a special 
problem of intermediate motions has been considered numerically. 


2. Basic equations 

Let E, H, and j denote the electric and magnetic field-strengths and the 
current-density, all measured in e.m.u. In the two-dimensional case, E and 
H are derivable from a vector-potential A, parallel to Oz, by the relationst 


H = curlA, E = —éA/ét. (2) 
Also if « is the electrical conductivity, 
curl H Aj 4ok(E+v AH). (3) 
Substituting (2) into (3), the magnitude of A is found to satisfy the equa- 
tion ' 
™ Ke A | 
V2A 4k _(v.grad)A}. (4) 
| et 


This is the equation valid inside the conducting cylinder; outside it 4 


satisfies vV24 — 0. (5) 


On the surface of the cylinder, A and its derivatives have to be continuous. 


+ In three dimensions it is necessary to include a term —grad V in the relation giving E 
In the two-dimensional case, however, E is a vector parallel to Oz but independent of 2; 
thus grad V can only be a vector parallel to Oz, independent of all the spatial coordinates. 
The contribution from grad V can accordingly be included in the —éA/ét term without 
affecting the relation div A = 0. 
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The steady motion v in a uniform liquid can be expressed in terms of a 


stream function % by the equation 


Vv ceurl(kys) = k A gradu, (6) 
where k is a unit vector parallel to Oz. Also in a normal mode of decay, 
Aa e~”; thus equation (4) becomes 

V2A tox{—oA-+k.(grady% A grad A)}. (7) 


This equation, together with equation (5) valid outside the cylinder and 
the necessary boundary conditions, determines the eigenvalues o and the 


eigenfunctions A. 


3, Small velocities 


When v is small, perturbation methods can be employed. When v = 0, 


equation (7) becomes 


V2*A+4aKoA 0. (S) 


Functions A satisfying this equation inside the cylinder, equation (5) out- 
side, and the appropriate boundary conditions, exist only for positive real 
eigenvalues o;. Let A,, A, be eigenfunctions (which may be assumed to be 
real) corresponding to different eigenvalues o,, ¢,. Then A, and A, can be 


shown to satisfy the orthogonality condition 


r s 


( ( A.A.dS 0, (9) 


the integral being over the cross-section S of the conductor. If there is 
degeneracy, the eigenfunctions corresponding to a given value of o can be 
chosen such as also to satisfy equation (9). Eigenfunctions are also nor- 


malized to satisfy 


[{ At? dS = 1, (10) 


r 


In solving equations (7) when ¢ is small but non-zero, suppose first that 
the unperturbed eigenvalues o, and eigenfunctions A, satisfying equation 
s) are all distinct. The corresponding perturbed quantities are then 


f xpressed iS series 


{ ite > a,,A,te > b,,A,+..., 
: (11) 
Co Ep, t e*, 


where ¢ isa small constant to which ys is assumed proportional. Substituting 


Ll) into (7), the terms of orders « and e? give 


> a,.(¢,—0,)A,+p,A, e1k.(grad¢% A grad A,), (12) 


s 


> 5,.(¢,—0,)A,+p, > a,,A,+v,A, = e1K.(grad¢% A grad > a,,A,). (13) 
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Equation (12) gives the first-order perturbation. Multiply this equation 
by A, and integrate over S, using (9) and (10). Then 


m, = fe ( [ k.(grad¢ A grad A?) dS 
se? | | K.curl(A? grad) dS 


de J A? grads. ds (14) 
5 
using Stokes’s theorem, where the last integral is a line-integral round the 
boundary C of S. But C is a streamline % = constant of the motion: thus 
the integral vanishes and p, = 0. That is, to the first order the motion 
does not affect the decay constant o. 


To get a,, (r 4s), multiply (12) by A, and integrate over S; then 


,.(F,.— Gg) e1 || A,k.(grad¢ A grad A,) dS. (15) 
The coefficients a,, are symmetric in the suffixes r and s; for 


rs s? 7 s 


(a a...)(o,—o,) « | | k.(grad¢ A grad A, A.) dS 


is A,A,gradys.ds 
C 
which vanishes like (14). 
The second-order perturbation equation (13) can be solved similarly to 
give v, and b... 


Using the symmetry property a,, = a,,, the results can be 
put in the form 2 
O,)O ings (16 


b..(o,. G.) VS (co O,)Ajy Aye (17) 
t 


Suppose that o, corresponds to the mode of slowest decay, o, to the next 
slowest, etc.; then o, < o, for all s, and so v, > 0. Thus the rate of decay 
in the lowest mode is accelerated by the fluid motion, to the second approxi- 
mation. No similar conclusion is possible for the higher modes, for which 
the series in equation (16) contains negative terms. However, no general 
reduction in the rate of decay is possible, since from equation (16) and the 
symmetry condition in a,, it follows that the quantities 


are all positive. 
The case of simple degeneracy can also be discussed without difficulty. 
If two real eigenfunctions A,, A’. correspond to the same eigenvalue o,, the 


zero-order eigenfunctions for the purpose of perturbation theory must be 
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taken as A.+7A’. There is now a first-order perturbation in the eigen- 
values o, given by 
i || A,k.(grady A grad A;) dS. (18) 


Ep 


’ 
Since this is purely imaginary, the rate of decay is still unchanged to the 
first order, but now the decaying field has something of the character of a 
progressive wave. To the second order, the effect of the perturbation is to 
nerease voth values of o by the same real quantity «*v,, where again r, is 


positive for the slowest mode of decay. 


4, Applications of perturbation method 


A simple application is to the problem of decay in an infinite slab between 


parallel plane walls y a. Here the unperturbed eigenfunctions and 
eigenvalues are given by ; 
. ° cos 
A (C cos py+Dsin py) . (ax), (19) 
sin 
4nKko p* |? (20) 
where either 
‘nt , D 0, tan( pa) x/p | 
; (21) 
( (0). cot( pa) x/ p| 


There are thus two interlacing sets of eigenvalues for any «a; the larger ones 


ire given approximately by 


pa vi7r/2, 4rKo x2 par? /402, (22) 
where v is a sufficiently large integer. Corresponding to the degeneracy 


indicated by equation (19), the zero-order eigenfunctions in perturbation 
theory must be taken as 

A = f(y)e*™, (23) 
where f(y C' cos py+-D sin py. 

Consider first the case u& ed(y), when the motion is entirely in the 
r-direction (parallel to the faces of the slab). To avoid infinities, the 
integrals in equations (9), (10), ete., are limited to a region S of the slab 
whose length in the x-direction is the wavelength 27/« of the eigenfunction 
in that direction. The first-order perturbation in oa, is found to be 


Ep, T taV, 


a 


| vf? dy 
where V “a 9 (24) 


a 


| f2 dy 


a 


«dd/dy) being the velocity of streaming. Thus the dependence of 
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the field on # and't is through a factor e+***-"®, This shows that the field 
drifts in the x-direction with the velocity V, which is a weighted mean of 
the fluid velocity. 
Secondly, consider a motion periodic in x such that 

ys = ed(y)cos(Bx+n) (25) 
with d(a) = 0 = ¢(—a). For simplicity, « and 8 are assumed to be com- 
mensurable; the integrations in equations (9), (10), ete., are taken over a 
length in the x-direction which is a multiple of the wavelengths 27/«, 27/8. 
One then finds that, in spite of the degeneracy, the first-order perturbation 
in o, is always zero. Thus a succession of oppositely-directed circulations 
does not, to the first approximation, lead to any drift of the field in a 
normal mode. 

The results for decay in a circular cylinder of radius a are similar. In 

terms of polar coordinates r,@ the unperturbed eigenfunctions are 

A CJ, (Kr) °° (n8), (26) 

sin 
where 7 is an integer, and K (= ,(47«o)) is a root of the equation 
J,_,(Ka) 0. 

In non-uniform rotation with angular velocity w(7) at a distance r from 
the axis, the first-order perturbation indicates a field rotating with angular 
velocity Q, where Y 
| ro(r)| J, (Kr) ]? dr 


Q = * _ ; ——» 27 


| r[ J, (Kr)? dr 
0 
so that Q is a weighted mean of w(r). On the other hand, a succession of 
oppositely directed circulations given by 
os ed(r)cos(mé+-») (m ~ 0) (28) 
produces no first-order perturbation of the eigenvalues, i.e. no migration 
of the decaying field. 


5. Fast motions: an exceptional case 

In the presence of a fast shearing motion, two kinds of normal decay are 
in general possible. In the first, there is little interaction between the 
motion and the field, the lines of force being almost identical with the 
streamlines. The decay tends to create a pattern of lines of force differing 
slightly from that of the streamlines; the effect of the motion is simply to 
push them back and so restore the original form of the field. Thus the rate 
of decay is similar in order of magnitude to that in the absence of the motion. 
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In the second kind of decay-field the lines of force run across the stream- 
ines: the effect of the motion is to elongate the lines of force in the direction 
fthe motion and so to bring elements containing oppositely directed fields 
nto close proximity. Thus in this case the motion greatly increases the 
ite of de ay of the field. 

The first type of decay . which includes only a limited set of modes, will 


considered in this section. In modes of this type, A is given by 


A A,()+-A,, (29) 
where A,(ys) represents the main part of the field, with lines of force 


coinciding with the streamlines, and A, is a small correcting term. If this 
expression is substituted in equation (7), the part involving A, disappears 
in the vector-product term on the right; in the remaining terms A, can be 


neglected compared with A,. Thus one finds that 


V?A, i7x{—oA,+k.(grad¢ / grad A,)}. (30) 


Since the vector product depends only on the component of grad A, along 
the streamlines, if ds and dn denote elements of length along and normal 


to a streamline, this equation can be written in the form 


ods OA 
V2A,+42KoA, Amn 1 (31) 
On Os 


Integrate this equation over the area 5S bounded by the two adjacent 


streamlines us, &+-dus. Then if C(ys) denotes the streamline ys, the right-hand 


side of equation (31) gives 
| aA, 
tirK OU ds dark dus | dA, 0. 
( ( me Cth) 
\lso, by Green’s theorem, 
| | v24, as | Oe gs | Othe ig, 
on J on 
j S ( , or Cd) 
d r OL ans L 
Sal . 10 Ag us : (Pw) Ay 
dis\ | On dis dfs , 
( (ub) 
where P(us) | OY ds. (32) 
J on 
Cd) 
Again || A,dS Ay Q(x) du, 
eran 1 
where (us) | (| ds, (33) 
con 


Cid) 
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Thus the result of integrating (31) is the (self-adjoint) ordinary differential 


equation 


l dA, 
Pata ib -4arKko Q(y)Ay 0. (34) 


This equation, together with the boundary conditions 
A, = 0, dA,/db = 0 


on the surface of the conductor, determines the eigenvalues o. Numerical 
values of o are most simply obtained from a variational (Rayleigh) prin- 
ciple; they are the stationary values of the Rayleigh quotient 


: (47k) | PUMy| it) aH] | @ Q(b) Az dib (35) 


as the function A,(s) is varied. An alternative expression to (35) is 


Wd 4.\2 a 
py (47x«)-} | | ( °} oradys 2 dS | | A« ds, (36) 
Jee diss , Pe 


the integrations again being over the cross-section of the cylinder. The 
advantage of (36) is that it frequently can be applied when explicit expres- 
sions for P(s), Q(s) cannot be found. 

Numerical results have been obtained for motion in the circular cylinder 
r = a, with the stream function 4 = Ky, where 


bo r(a2—r?) cos 6/a8. (37) 


This represents a double circulation of the fluid; to give rapid motion, the 
constant K must be large. Equation (36) is used, with the approximating 
function A, vh2+-Bus4. The least value of © then gives ozxa? = 4-41 as 
an approximation (probably a good one) to the lowest eigenvalue. For 
comparison, in the absence of motion the lowest eigenvalues corresponding 
to n l and n 2 in equation (26) make oz«a? equal to 1-45 and 3-67 
respectively. Thus the decay is speeded up by the motion, but remains of 
the same order of magnitude. 


Fast motions: the slab problem 

The second type of decay in the presence of fast motions (that in which 
the lines of force run across the streamlines) will be studied in two special 
cases. These are non-uniform streaming in an infinite slab, and non- 
uniform rotation in a circular cylinder. 

Consider first the slab problem. As in section 4, the slab is bounded 
by the walls y 


La, the motion v between the walls being in the 2- 
direction. For simplicity we consider a linear velocity-profile, given by 


v=py, w= —fuy?. (38) 
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Then equation (7) reduces to 


V?A+4rKoA 4rKkpy = : (39) 
\ solution periodic in x is taken to be 
A= fig“, (40) 
where f satisfies the equation 
d} ° . ; 
he (4arK0—a*—4rKpray)f = 9. (41) 
\t the walls of the slab, A must join up smoothly with solutions of 
Laplace’s equation which vanish at y too. These are 
A Cexiz—-uv)-t (y >a), 


Deriz+u-ot (y <a). 
Thus, to make A and its derivatives continuous at the slab walls, f(y) must 
satisfy the boundary conditions 
f(a) vf (a), f'(—a) xf(—a). (42) 
Non-dimensional combinations of the different parameters are defined 


DY the equ tions 


q = aa, p = 4nxaryq | (43) 
w p-\(4rKa*0—q’), Y y/a 
Then equations (41) and (42) become 
df sie 
7y2" p(w—tY)f 0 (44) 
nd f’() gf), f(—l)=af(—)), (45) 


where in equation (45) f is regarded as a function of Y. The fact that the 
motion is fast implies that p is large; the magnitude of w depends on 
that of o. 

Let new variables Z, X be defined by 


= p(w—iY), X 50Z}/p. (46) 

Then the solution of equation (44) is expressible in terms of Hankel functions 
in the form z 4 a ” 
— f — Z\[ PH®(X)+QH2(X)]. (47) 


On inserting this expression into equations (45), and using the recurrence 
relations for Hankel functions, one finds that 


Z\[PH%(X,)+ QHE(X,)] = —q PHY(X)+QHP(XD] gg) 
Z} PH 1) (X4) QOH : (X,)| q| PH\)(X,4) T QH(X3)| 
where r Z 
Z,=pw—i), Z,=plw+i) | (49) 


X, si p?(w—2)?, X, = Rip'(w- iyi) 
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If the ratio P/Q is eliminated between equations (48), the result is 
(Z, Z,)*|H™.(X,)H®;(X,)—H@,(X,)H™(X,)] 
Zi go HY,(X HY (X,)—H(X,)H(X,)]4 
+ Z3 q[ H\?(X,)H®,(X,)—H?(X,)H™,(X,) | 
q*| HS? (X,)H\(X,)—H?(X,)H{P(X5)] = 0. (50) 


This is an equation giving the eigenvalues of w, and so of the decay con- 
stant o, for given values of p and q. 


7. Asymptotic expansions 

By equations (49), Z, and Z, are in general large of order p, and X, and 
X, of order p?. Thus in equation (50) the Hankel functions can be replaced 
by their asymptotic expansions. Since for all values of X these indicate 
that H,(X), H®.(X) are comparable in order of magnitude with H{)(X), 
H‘*)(X), the term involving (Z, Z,)! in equation (50) dominates the other 
terms, so that the equation can be simplified to 


H”,(X,)H@.(X,)—H®,(X,)H%,(X,) = 0. (51) 


The leading terms of the asymptotic expansions are as follows, for 
different ranges of arg X: 


9 \1 
H\))(X) | ze X—tu7-t7)  (—7 arg X 277) (a) 
2 \1 ~ ; : 
_ ( a {2 cos v7 (X+v7+}7)_ pi(X —-ly7- i) ( In < arg xX < 0) 
(b) 
9 \1 
H\?)(X) ( a (X-jvnr-it) (—Inr < arg X < 7) (c) 
= y *£2 cos vz eX thya+im)_| p—i(X—hy7 im) (0 < arg a 4 27). 
(d) 
In the range —7 < arg X < 0, where both (a) and (b) are valid, the extra 


term in (b) is negligible; a similar remark applies to (c) and (d) in the 
range 0 < arg X < 7. 

Using these and similar relations, equation (51) can readily be shown to 
possess solutions only if w has a positive real part, and an imaginary part 
lying between —1 and +1. When w lies in this range, 


- . - r 5 mt) 
—\|n < arg X, < 4m < argX, < 5z. (52) 


The expansion (a) can accordingly be used for both H(X,) and H{?(X,), 
i a v ( 1 L 
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i 


(c) must be used for H‘*)(X,) and (d) for H'?(X,). On use of these 
xspressions, equation (51) becomes, after simplification, 

exp{i p!(w+-7)!!—expf{s p?(w—7)!} . (53) 
fhe solutions of this equation fall into two groups, corresponding respec- 


vely to real and complex values of w. 


§, Real eigenvalues 


If w is real, let us write 


(ao 1 )? pe 10 . (54) 
here y and 6 are real. Then equation (53) further simplifies to give 
sin(4p?d) 5 exp( 4 py). (55) 


In this equation, y cannot be appreciably negative, since this would make 
the right-hand side large, whereas the left-hand side cannot exceed unity. 
Now y = 0 corresponds to w = wy = 1/v3; hence the possible values of w 
re effectively limited to the range w > w». By equation (49), this in turn 
mplies that 47«xa*o is large of order p, so that decay of the field in these 
modes is fast compared with that in the fundamental mode in the absence 
f motion. 

Ifw is appreciably greater than wy», p'y is large, and the right-hand side 
fequation (55) is negligible. Thus in this case approximately 

sin({p'd) = 0, = 4p'd = vz, (56) 
here vy is a sufficiently large integer. In particular, if the decay is very 
pid, w is large and 8 su, so that 


2(wp) VT, 4oKko x“ by? 2, (57) 
Equation (57) gives the same expression for the decay constant o as 
equation (22), valid in the absence of fluid motion. This is not surprising; 
{the decay is sufficiently fast, the motion has insufficient time to modify 
t appreciably. 
All save one or two values of w are given very closely by equation (56). 
But in this equation 5 cannot be less than (4)!, its value for w = wy; thus 
vr > (4)ip?. (58) 
this implies that as p increases the number of modes of decay with a real 
continually decreases. In fact, the number of intersections in the (x, w)- 
plane of the two curves 
x sin(4 70), x Fexp( 4 pty) 
continually decreases as the increase in p compresses the sine-curve towards 
the axis w — 0, and carries progressively more and more of its arches past 
the value w — w, for which the exponential is unity. One can expect that 
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the eigenvalues of w will coalesce in pairs and then be replaced by pairs of 
conjugate complex eigenvalues. This conclusion is confirmed by the 
numerical work described later. 


9. Complex eigenvalues 


i] 


Consider now the complex solutions of equation (53). The assumption 
(52) made in deriving this equation corresponds to 
ba < arg(w—i) < 0 < arg(w+i) < 3 


In Fig. 1 the lines AB, CD are the lines 


3 


l 
7 
3 


arg(w—?) —in, arg(w-t7) 
in the complex w-plane; on A B, (w—7)! is purely imaginary; on CD, (w+i) 


is purely imaginary. If wis appreciably above A B, the term exp{4p!(w—i) 


i,A D 





O 








-t"C B 


Fic. 1. Lines in the complex w-plane. wy is the point w 1/v3 


of equation (53) is large. Thus in this case equation (53) can be satisfied 
only if w is appreciably below CD, so that expf{4p!(w--7)!} is also large. 
But in this case the equation becomes 
exp{ip!(w+i)!} = exp{4p!(w—i)}, 

which is satisfied only by real values of w in the ranges of arg(w--7) under 
consideration. Thus no complex roots of equation (53) exist such that w 
is appreciably above AB; similarly, none exist such that w is appreciably 
below CD. 

If w is appreciably below AB and above CD, both of the exponentials 
on the left of equation (53) are small, and their difference cannot equal /. 
Thus the only possibility for complex solutions of (53) is that the point v 
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should be very near one or both of the lines AB, CD. The intersection E 
of the lines is the point w = w, near which the transition from real to 
complex eigenvalues takes place. Only a small number of solutions of 
53) can be expected to correspond to points w near to E; most of the 
omplex eigenvalues correspond to points near AE or CE, but well away 
irom E. 

Consider, for example, values of w corresponding to a point near AE 
but appreciably above CE. Then the first exponential on the right of 


equation (53) is negligible, and the equation becomes 


exp} 4 p?(w—z1)*} t. 
‘ g \l , . 
This has roots Ww d | } [ev t } )ar |i st (59) 
L6p Ps 


where v’ is a positive integer or zero. Similarly, the roots of equation 
53) corresponding to points near CE but appreciably below AE are given 


; g \i ae 
uw 2 | }? [(2» } 1 Jar |e Te (60) 
L6p 
The condition that the points w given by equations (59) and (60) must 
it reach F is that 7 4 
(21 bar < (§)‘p?. 
Thus the total number of complex eigenvalues is roughly 7-1(4)‘p!. This 
grees, as is to be expected, with the number excluded from the set of real 
eigenvalues by equation (58); effectively the real eigenvalues which have 
lisappeared have been replaced by an equal number of complex ones. 
The lowest eigenvalues given by (59) make w—i small, of order p 
This means that the asymptotic formulae used in deriving equation (53) are 
not strictly applicable to the lowest one or two of the eigenvalues; in fact, 
59) gives , ae p ’ 
, X, s0p?(w—?)? }(2v + t)z, 
which is not large for the lowest eigenvalues. However, equation (58) should 
give all save the lowest one or two eigenvalues satisfactorily, and even for 


these it should at least give the correct order of magnitude. 


Let w = c-+-id, where c and d are real. Then, by equations (43), 
oC re rA( C4 id) t x2 (47x). (61) 


The decay field is proportional to e’®*-%, and so includes a factor 
pro} 


eae padi) 


This implies that the field is undergoing a translation with velocity pad. 
For the eigenvalues given by equation (59), 0 << d < 1; thus the velocity 
is that of the fluid at a point y = ad between the central plane y = 0 and 
the boundary y = a. Similarly, for the eigenvalues given by (60) the field 

















398 T. G. COWLING AND A. HARE 


moves with the velocity of a point y = —ad of the fluid below the plane 
y 0. For the lowest eigenvalues (v’ small in (59) or (60)) the velocity of 
the field is very nearly that of the fluid at one boundary. 
The rate of decay of the field is given by the real part of o, which js 
found to be 
(477) -1(a2-+ pe/a2) (darn) 2 »|'l. (62) 
| 2a*| 4 ‘ 
Thus the rate of decay increases with p, roughly like p'; for all the lower 
modes it is much faster than that in the absence of motion, given by 
equation (22). 
The transition from real to complex eigenvalues takes place at 


w Wo 1L/v3. 


This corresponds to a time of decay o~! roughly equal to 

47Ka*\3/p V3/(uad). 
Since the fluid at the walls y La travels a wavelength 27/« of the field 
in time 27/ (yaa), real eigenvalues are impossible if fluid at the walls travels 
more than about a quarter of a wavelength in the time of decay. When 
the motion is faster than this, the field is carried about with the fluid. 


10. The decay fields 

In modes corresponding to real values of w, the decay fields can be 
expected, like the eigenvalues, to be little affected by the motion. The 
decay fields for the complex eigenvalues have more interest. 

Consider values of w corresponding in Fig. 1 to points of AE not close 
to EF. The field corresponding to the vector-potential of equations (40) and 
(47) is 


H, —(Z/a)[ PH™,(X)-4 QH®,(X)]ei-ot, 
H, inZ*| PH\)(X) | QH\?(X) |e iax—ot 


Here the ratio P/Q is found to be e!'". Hence, using the asymptotic forms 
(a) and (ec) of section (7), one finds that, except for values of y very nearly 
equal to a, 


H, QZ ( x) i(X+7) piax—ot 
oe (63) 
H, QuZ || x)" i(X+%#7) piax—ot 
7X 
In this, X 3ip'(w—iy/a)!. Because of the factor p!, the imaginary part 


of X varies rapidly as y varies, and it has a maximum value when 
y/a im(w) d. 


Thus the magnetic field attains its maximum value on the plane y = 44, 
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where the fluid and the field have the same velocity. Well away from 
this plane, the field is negligibly small. Thus it is always small on the face 
a of the slab; it is also small on the face y = a, save for the lowest 
few eigenvalues, for which the maximum field is attained near y = a. 
From equation (63) the ratio H,/H, is of order Z, and so is large of 
der p! (save near the boundaries of the slab). This shows that the loops 
nthe lines of force are greatly elongated in the x-direction as a consequence 
f the motion. 
The case when w corresponds to a point of CE in Fig. 1 can be treated 
similarly; the only difference is that the maximum field is attained below 


the central plane y 0. 


11. Other velocity profiles 

In the problem just considered, the existence of a family of real eigen- 
values depended on the skew-symmetry of the velocity-profile with respect 
to the central plane. No such family would exist for a less special profile, 
though the eigenvalues for a general velocity varying monotonically across 
the section can be expected to have properties generally like those for the 


) 


ar ve locity profile. 


The properties of eigenvalues when the variation of the velocity is not 
monotonic can be roughly inferred from those for the ‘roof-top’ profile 
iven by 

(a—y) (0 y a), v plat+y) (-—axy<09). (64) 
The discussion of this profile exactly follows that for the linear profile; 
ience only the main results need be quoted. In the sufficiently high modes, 
the field is carried along with the velocity of the fluid at y + 4a; the 
eigenvalues for very rapid decay are given by 
iaKo 72+ dy*a-2+ (4a) duaar, (65) 
which differs from equation (57) only through the imaginary part which 
represents the velocity of translation of the field. 

The lower modes (those corresponding to the complex eigenvalues of 

section 9) fall into two groups. In one group, the field moves with a 


velocity less than }ua, the velocity becoming small for the lowest modes. 


In the other group, the velocity of the field is greater than }ya, and 
ipproaches ua for the lowest modes. The field in either case attains a 
strong maximum on the planes where it moves with the velocity of the 
lluid; these planes are near the walls of the slab for the first group, near 
the central plane y = 0 for the second group. The first group corresponds 
to pairs of virtually) equal eigenvalues, given by 


torKo 2+ 3(1+-tv3){3(2v4+-$)2? pax}. (66) 
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The appearance of these pairs of values corresponds to the fact that the 
field is vanishingly small near y = 0, so that effectively decay modes exist 
such that the field is limited to one or the other of the regions y > 0, y < 0, 
No such pairs exist for the second group, for which the eigenvalues are given 


by 4aKo v4 4p rai +43(1—iv3){3(v+43)z2uaK}4. (67) 


As regards other velocity-profiles, when the velocity is an analytic 
function of y, general methods of the theory of functions of complex 
variables can be used to determine the eigenvalues (4). Such methods are. 


however, of limited application here, and will not be discussed. 


12. The cylinder problem 

Non-uniform rotation in a circular cylinder can be discussed like the 
streaming in a slab. For illustration, consider the case when the angulai 
velocity at a distance r from the axis of the cylinder is Ar. In terms of 
polar coordinates 7, 6 equation (4) then becomes 


4 4 1cA ' c2A ane“ | wr), (68 


cr r or ro? Ot 00, 


Assume a solution of form 
A r firye ind ot 
where » is a positive integer. Then the radial function f(r) must satisfy 
d?f ; n2—1) , 
“ot | deracto nar) . ay 0. (69) 
dr? * | ye | 
This has to be solved subject to the boundary conditions that f is finite 
at r 0, and df/dr (n—4)f/r at the boundary r = a of the cylinder. 
[f A is large, we may expect from the results of the slab problem that ¢ 
is in general of order nAa. Thus, within a range 0 < r < (4z«Ad)!, the term 
involving mAr in equation (69) is negligible compared with the o-term; in 
the range (47xA)! < r < a, the term (n?—})/r? is negligible compared with 
the c-term. Approximating by neglecting these terms in the appropriate 
ranges, we can obtain asymptotic expressions for the solutions of equation 
(69) in the two ranges; both asymptotic expressions are valid near 
r (47xA)-*, 
which enables the arbitrary constants in the two to be linked up. 
When this is done one finds that, in order to satisfy the boundary con- 
ditions, an equation must be satisfied of the form 


exp{$p? (w+)! —ina—tir}—expf{f p'(w—z)!} i, (70) 
where p = 4nxa'nd, w = —i+2e/dan. (71) 


Kquation (70) corresponds to equation (53) for the slab problem. 
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\s in the slab problem, there are two sets of solutions. For the modes 
f very fast decay, w is real; this implies that the field rotates with the 


ngular velocity 3Aa of the fluid at r = 3a. For modes of less rapid decay, 


the field rotates faster or slower than this; for the low modes, it rotates 
with the angular velocity of fluid near the axis of the cylinder or near the 
surface. The decay field is small save where the field and fluid are rotating 


with the same angular velocity. 


13. Moderate velocities 
In order to indicate how the results for low and high velocities join up, 


he slab problem of section 6 has been studied numerically for moderate 


ocities. The range of velocities considered is that in which takes place 
he transition of the lowest pair of eigenvalues from real to complex values 
The equation to be studied is (44), i.e. 
d*} -»Y)f 
(p*—ipy ) 0, (72) 
jy | pt )f 
here p? — pw; the boundary conditions are 
f’(1) qf (1), f'(-—)) = @af(—)). (73) 
[f p? is real, any function / satisfying these conditions can be shown to be 
I the tol . ¥ on oo 
f = e4(U+%iDV), (74) 
ere U’ is even and V odd in Y, and ¢ is a constant. This implies that, 
Lan ire the real and imaginary parts of f, 
0)/v'(0) e(O)/u(O), (75) 


h of these expressions being in fact equal to —tan¢. 


\ccordingly the calculation for real eigenvalues proceeded as follows. 
For simplicity the quantity q¢ xa was taken to be unity. A solution satis 
tying the boundary condition f’(—1) f(—1) was begun at Y 1, the 
Taylor series at Y | being used to give the values at Y 0-9 and 

0-8. The values at } 0-7, —0-6, —0-5...., were then found by 

peated use of the approximate formula 

hf "la l(a h) f(a h) 2f(a) (76) 
h / 0-1. This formula neglects terms of order /*, but is sufficiently 


urate for present purposes; its use is far less laborious than using the 
Taylor series. Corresponding to a suitable value of p, solutions were 

bulated for several different values of p?. The correct value of p? is that 
lor which the solution satisfies equation (75). 

The variation of p? with p is shown in Table 1 for the lowest two modes. 
The values corresponding to p 0 are those given by equation (21). As 
pd 
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p increases, the two values of p? come together, and finally coalesce at 
p = 3-6, which therefore corresponds to the transition to complex eigen- 
values. This value of p implies that fluid at the slab walls travels about 


one-sixth of the wavelength 27/« of the field in the decay time o~!. 


TABLE | 


Complex eigenvalues are less easy to discuss, because of the absence of 
a symmetry condition like that of equation (74). The method finally 


adopted involved abandoning assigning a value to q beforehand. Put 


p? = q(y+2), Y = q-'y+5/p, QO = pi¢’. (77 
Then the equation satisfied by f is 
d*f dy } (y—tiQy)f 0 (78) 
and the boundary conditions are 
f'(n)/f™) I at 7 q(1+6/p) = 6, say |) (79 
—l at 7 q(1—9d/p) c, say ” . 


If f, and f, are two independent solutions of equation (78), the general 
solution is f = Af,+ Bf,. On elimination of the ratio A/B between the two 
equations (79), one gets 
SiO)+f6) — file)—flo) 
flb)+fo(b) — f3(e)—fele)’ 


Accordingly, to obtain complex eigenvalues assigned values were adopted 


(SU 


for Qand y. The two independent solutions /,, f, were obtained by numerical 
integration of equation (78). The (complex) functions on either side of 
equation (80) were then tabulated for negative values of b and positive 
values of c. If Q was large enough, values of b and c could be found graphi- 
cally to satisfy equation (80); knowing these and Q, values of q, p, and 4 
could then be derived. 


When q = 1, the transition from real eigenvalues to complex was found 
to occur at Q = p = 3-6, p? = y = 2-7. Thus to get complex eigenvalues 
the values Q = 4-0, y = 2-7 were adopted; this gave 

q 0-99, ra) Q) +0-16. 


The value of 6 implies that the field moves with an x-velocity equal to that 
of the fluid at y -0-16a. Since q differs so little from unity, the results 
so obtained can be regarded as continuing those for real eigenvalues with 
q b. 
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Lines of force for both real and complex eigenvalues are shown in Figs, 
2-9, where the curves correspond to different constant values of the real 
part of A. In each diagram, the upper and lower boundaries are the 
walls y +a of the fluid; the breadth is one-half of a wavelength of 
the field. It can be seen that increased values of the velocity parameter p 





\ 











Fic. 8. Lines of force corresponding to p 3-6, 


for which first and second decay modes coalesce 











N 





Fic. 9. Lines of foree corresponding to Q 4, 
2-7, for which fields in the first pair of decay 
modes are being transported with the fluid 
at first lead to a stretching of the lines of force in the direction of motion, 
and to a closer approach of the patterns representing the first and second 
modes to each other. After the coalescence of these two modes, complex 
eigenvalues appear which correspond to a good deal of asymmetry of the 
field with respect to the plane y = 0. 


14. Conclusions 
The results outlined above make it possible to predict the nature of decay 


fields in the presence of an arbitrary two-dimensional motion. Such 4 
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motion can be divided up into circulations in separate regions, bounded 
by streamlines passing through stagnation points. In the interior of any 
such region, the decay fields will vary roughly as in the problems con- 
sidered above. In slow motion, or in fast motion when there is very rapid 
decay, the field in any region will be carried round with a period of circula- 
tion which is a mean of periods of circulation in all the regions. 

In fast motion when the decay is less rapid, the difficulty of transporting 
es of force through the stagnation points can be overcome only by 
supposing that very few of the lines of force cross the boundaries of the 
regions, and that these are mainly concentrated just in front of the stag- 
nation points. In this case the lines of force, when not nearly coincident 
with the streamlines, will be greatly extended in the direction of the motion. 
The field will be transported by the circulations, and will attain maximum 
strength on streamlines where its period of circulation equals that of the 


fluid. The modes of slowest decay will circulate with nearly the period of 
circulation characteristic of streamlines on which this period is a maximum 
a minimum. These may include streamlines immediately surrounding, 


r passing through, a stagnation point. In the first case the decay modes 
n effectively be taken as modes in which the field is confined to one of the 
the second, when the field circulates with the fluid near the 


boundary of the regions, some slight interdependence of the regions is 

ssible. As the velocity increases, the period of decay in the lowest modes 

ries roughly inversely as the 2 power of the velocity (save for the special 
des for which the lines of force and streamlines nearly coincide). 
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NON-UNIFORM SHEAR FLOW PAST CYLINDERS 
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SUMMARY 

A general method is described whereby an approximation of any desired degree 
of accuracy to the stream functions for two types of variable shear flows past 
finite cylinders can be obtained. The two shear distributions in the free stream cai 
be approximated to the linear shear distribution and the shear present in an 
unretarded incompressible boundary layer respectively. In every case the stagnation 
streamline is displaced from the position opposite the line of symmetry of the 
cylinder, and general expressions are obtained for this displacement. The line of 
symmetry may be in the direction of or perpendicular to the direction of flow. 
The two particular examples cited are those of a general elliptic cylinder and 
cylinders of the form r = p(1+y cos 8), where |y| < 1, r and @ being the polar 


coordinates, and 2p the maximum width of the cylinder. 


Nomenclature 

al Bessel function of the first kind. 

I Modified Bessel function of the first kind. 

Y Bessel function of the second kind. 

K Modified Bessel function of the second kind. 

l Typical linear dimension. 

L Perimeter of the cylinder. 

N Rotational flow parameter (—= U/wo/l). 

p Typical length in the cylinder. 

q Velocity. 

q Velocity on the cylinder. 

G.: Cy Velocity components in the x#- and y-directions. 
dr» V0 Radial and transverse velocity components. 

r, 0 Polar coordinates. 

U Standard velocity. 

x,y Cartesian coordinates. 

Y. Deflexion of the stagnation streamline in the free stream. 
x Angle the tangent to the cylinder makes with the x-axis. 
y General parameter. 

‘a Circulation. 

5 Thickness of the incompressible boundary layer. 
€ Eccentricity of the ellipse. 

us Stream function. 


[Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 4 (1957)] 
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Disturbance stream function. 
Vorticity. 


Laplacian operator. 


1. Introduction 


Tue stream function & for the steady, rotational two-dimensional flow of 


n incompressible inviscid fluid satisfies 


2), 2d 
CW Cc~W 
-+—+w 0, (1) 
os cy~ 
where the vorticity w, given by 
oC Cd... 
- ee 9 (2) 
Cx cy 
; : CW CW ‘ 
satisfies the equation q rq, 0, (3) 
Cx cy 
nd q <* ~ are the velocity components parallel to the z- 
7 Ca 


d y-axis respectively. 

In a previous paper, Mitchell and the author (1) obtained the stream 
functions for some variable shear flows past circular cylinders, and com- 
ared certain aspects of these flows with experimental results. As far as 
he author is aware, no theoretical solutions of (1) and (3) exist, with the 
xception of those obtained in (1) and those for constant shear flow. In 
the present paper, a method is derived for obtaining an approximation to the 
stream functions for some non-uniform shear flows past general cylinders, 

dthe method is outlined for obtaining the flows past an elliptic cylinder 


nd any cylinder whose equation is of the form r = p(1+-y cos 8). 


2. The variable shear distributions 
The general solution of (3) is 
Ww I (), 


ind so the vorticity is constant along a streamline. In the particular case 


whe 
— Ww f (ys) Cf, /]*, 


vhere / is some arbitrary length, equation (1) becomes 


Cs Ob “ 
oF in (4) 
ox™ cy 
The two cases incorporated in (4) give an approximation to a linear shear 
distribution taking the positive sign, and a shear distribution approximating 
to that present in an unretarded incompressible boundary layer, taking the 


hegative sion. 
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(a) The ‘linear’ shear distribution 
The appropriate equation (4) is 

Cus Cus 


> 


- os /1?, (5) 
Cx" cy 
A flow distribution parallel to the a-axis in the free stream near x = —zx. 
consistent with this equation is given by 
us aed! Abe yil (6) 
which on differentiating gives the velocity distribution as 
gd. I "ave yl ___be yl), (7) 
and the vorticity distribution as 
Ww l taev! be-v"), (8 
where the constants a and b are arbitrary. Introducing a standard rotation 
w, and a standard velocity U as the values of w and q, at y 0. it follows 
from equations (7) and (8) that 
a=WU—w,)), + 11(U +-w,)). 


With these values for a and b substituted in (7) and (8), the velocity and 


vorticity distributions in the free stream near x © become 

q,/U = eosh y/l—(1/N)sinh y/7, 9) 
and w/w, = cosh y/l—N sinh y//, (10) 
respectively, where V U/wol. In order to keep the velocity positive 
everywhere, the parameter V must satisfy the inequality V > 1. Expand 


ing the above distributions in terms of y//, and neglecting higher powers 
of y/l, they become 

q,/U = 1—(1/N)(y/I)+4(y/))? (11 
1— N(y/l) (12) 


and W/W» 


for the velocity and vorticity respectively, and so the free stream shear 
distribution (10) is approximately linear in the vicinity of the x-axis. 
(6) T he ‘boundary laye r’ shear distribution 

The appropriate equation (4) is now 


9 9 
\ c 2; Cus 


ox* Oy" 
A flow distribution parallel to the x-axis in the free stream near x —o 
consistent with (13) is 
ys = acosy/l+bsiny/I, (14 
which on differentiating gives 
q, = [bcos y/l—asin y/I], (15) 


and w = l-{acosy/l+bsiny/I], (16) 
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where a and 6 are arbitrary constants. If a standard rotation w, and a 


standard velocity U are introduced as the value of —w at y = 0, and the 


ue of g, at y = 5—y, respectively, it follows from (15) and (16) that 





a w, I, b = wyl*| N—sin(d—y,)/1|see(d—yp)/l, 
vhere .\ U'/wol, and yy is a variable parameter lying in the range 
y If 


COS Yy/l ] N, 
nd the arbitrary length / is chosen to be 


l (2/z)6 


he velocity and vorticity distributions in the free stream near a X 
ecome q [ sin 57(y -Yy)/O, (17) 
d a m N cos Sa(y Yo)/O; (18) 
respectively. Intherange 0 < y+y, < 4, the free stream velocity distribu- 
17) is Lamb’s approximation for the flow in an unretarded incompres- 

sible boundary layer of thickness 6 (2). 


3. The general cylinder with its line of symmetry in the direction 
of the stream 


\ cylinder is now introduced into the flow, in such a way that the line of 
metry lies on the x-axis, the origin lying within the cross-section of the 


ylinder, the equation of which may be written as 


? po(@), (19) 
here id @ are the polar coordinates and p is some typical length in the 
viinder. The function 4(@) must be an even function with no zeros or 

singularities in the rang 7 4 < 7. As the method is developed, it is 
seen that these restrictions are necessary. The problem, which is illustrated 
Fig. 1, is to find the stream function which satisfies in the first instance 
juation (5) and takes a constant value on the cylinder. Considering the 
inear’ distribution in the first case, and retaining the general a and 6 in 


ie function ’ is introduced by the transformation 
aevl it be-vl Lp, (20) 


From (5), 4’, which tends to zero as 2 tends to ©, satisfies the equation 


~ = P/P. 21) 
Cai” cy 
Using polar coordinates r and 6, and the method of separation of variables, 
21) is solved to give 


. 4 B, 6K, (r’) > (A, cosn6é+ B, sinn6)K,(r’), (22) 


n=0 
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where A,, and B, (n = 0, 1, 2,...) are constants, K,, is the modified Besse] 
function of the second kind, and r’ = r//. Thus, from (20), the required 
stream function % is given by 


bb = aeY!+be-Y'+ BOK, (r’)-4 (A, cosné+ B, sin n@)K, (7’). 


t 


n=0 














Ty 
||P 
| | : 
r= po(6) 


L=-O 
(Velocity profile) 
Fic. 1 


Now ys must take the constant value, say c, on the cylinder given by (19), 
and taking for convenience / — p, this gives 


c = aexp{d(@)sin 6}-+-b exp{—¢4(6)sin 6} + B, 6K, {d(0)} + 


> (A, cos né- B, sin nO) Kk, {d(0)}. (23) 
n 0 
The constants A, and B,, (n 0, 1, 2,...) are found by equating coefficients 
of sinné and cosné, the requisite functions exp{+4(@)sin 6} and K,,{d(8)} 
being expanded as Fourier series. Write 


exp{+4(@)sin 6} > X,,cosmO+ > Y,, sin mé, (24) 


m=0 m=1 


where X,, and ¥Y,, (m = 0, 1, 2,...) are evaluated by the usual methods. 
Appendix (qa) illustrates how this is carried out in the case of an elliptic 
cylinder, where r = p/(1+-ecos@), p being the semi-latus rectum and « 


the eccentricity. The Bessel function is written as 


K,,{6(8)} = A, 9 +A,,1 cos 6, (25) 


nt ’ 


where A,,, and A,, , are the first two terms of a series expansion of the form 


> A,,,cosy#, amended to give a good approximation to K,,{d(6)} in the 
0 : 


range —a7 < @ < 7, when this series is terminated at I. 
£ i 
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The expressions given by (24) and (25) are now substituted in equation 


93) resulting in 


(a+b) > X,, cosmé+(a—b) > Y,, sin m6+ By (Ag 9 +-Ao C08 8) 


0 m=1 


i 


> (A,, cos né+- B,, sin 6)(A,, 9 +A, cos 4). 


n U0 


[he constant ¢ in the above expression is as yet arbitrary, and in fact, as 


sseen from the subsequent analysis, remains so, unless use is made of the 
lowing hypothesis (as in the case described in (1)). The hypothesis is 
hat the introduction of a cylinder into an inviscid, incompressible, rota- 


l 


nal flow does not alter the circulation round a path L, which coincides 


with the perimeter of the cylinder. This assumption was used by Tsien (3) 


see also Richardson (4)) and Mitchell and the author (5) in the case of 
nstant shear flow past cylinders. Equating coefficients of the sine and 
sine terms, from which it is seen immediately that B, = 0, two sets of 
and B,. 


he equations, being of the following form, do not give direct values. Con- 


juations are obtained, whose solutions give expressions for A 


n 


sidering the cosine terms in the first instance, 


\n oA 11 A c—(a+b)X, 
No.1 A \1041+4A., 40 (a+-b)X, 
4A, A+ Aso Aot$A31 As (a+-b)X, 
4Nn—-1.14m—-1+AmoA Ans aA (a+b)X,, (m> 1). (26) 
\ similar series is obtained by equating the coefficients of the sine terms 
oon 
1,9 B,+4A,, Be —(a—b)Y, 
LA, Bi t+ Ao» B.+4A;3, Bs —(a—b)Y, 
tro, Bot Azo Bg+4A,, B, (a—b)Y, 
4A n—11 Bm-1tAmo 2B 1: Bass (a—b)Y,, (m>1). (27) 


he dependence of the solution on c is illustrated in (26), the A, being a 


lunction of the constant c. Using the hypothesis stated above, therefore, 


the constant ¢c is evaluated as follows. 


fhe stream function is given by 


yb = ae/P + be-v S (A, cosn6+ B, sin n6)K,(r/p), (28) 
0 


nm 
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from which the velocity components in polar coordinates are given by 


lL cus Cds 
r 00 OV 


dr 
The velocity at any point is given by g? = q?-+-qj. In particular, the velocity 
on the cylinder is given by q,, the velocity parallel to the tangent to the 
cylinder at its point of contact, as illustrated in Fig. 2. The velocity per- 
pendicular to the tangent at its point of contact is necessarily equal to zero, 
as may be easily checked. 


e 4 








The cylinder in cartesian coordinates is given by 
(a?+-y?)! pd(tan-ly/x), 

2 P dy A(@)cos 6-+-4'(8@)sin 8 
from which tan a ' 

dx 4(8)sin 0B—4'(0)cos 0 
with « as shown in Fig. 2 and 4'(@) d d6| (0) |. This relation defines 
sin « and cos a in terms of @ and 4(@). Resolving the velocity components 
parallel to the tangent at the point of contact 


t 7 ° \ » 

q 19, COS(a 7) +-qgsin(« 6)}|, pdl(6) (29 
with sina and cosa, obtained from above, substituted in equation (29). 
Stagnation points on the cylinder are given by q 0. The circulation I’, 


taken round the perimeter L of the cylinder is given by 


r=N+0,=4¢qadi, 


i 
where di is the element of length on the cylinder, and I, and [, are the 


circulations arising from the undisturbed flow and the disturbance flow 


respectively. The element dl is given by 
dl = [(dr)?+r°(d0)?}! 
pL{d’(4)}*+-{6(8)}*}! dd. 
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From equations (28) and (29), the circulation T, becomes 


2 r 1 cos( y 0) , 
; | S nA, sinnO+nB, cos n@)(A,, 9+A,, 1 cos 8) 
; aol) _— : 
x — { n 
sin(a—6) (A, cos n6-+- B, sin n@) (A, 9 +Ay.1 COS 8) + 
50 \A(@) : : 
1,0 +A, 140080) |[{4(0)}2+ (¢'(9)}?|§ dé, (30) 
sai ae alate 1K,,(r’) :, 
here in (30) A, {4(@)} and | have been replaced by the approxi- 
dr S(A) ‘ 
r @ 
\tions given by (25), since the velocity q, is that velocity on the cylinder, 
th sina and cosa replaced from the above expression for tan a in terms 
#. The condition that [T, = 0 from (30) gives a relation involving only 
ind c, where the A,, are replaced by the appropriate known functions of 
ind c from (26), resulting in an expression of the form 
( g(a, 6, X ns Amor Ami» Ao) (31) 
here g(A,) is a known function. 


[f this value for ¢ is substituted in equations (26), they become 


) 1, (A, oA, 1 A, (a b)X, 

Ao B Ai0 i. tr, Ao (a b)X, 

AX n—1.1Am—1 tAm.oAmt+3Am+11 441 —(a+b)X,,,. (32) 
he systems of equations (27) and (32) are of the same type. It is seen that 
f either system is terminated at a given value of n, it contains (n-+-1) 
nknowns and only x equations. Thus, unless the number of equations is 
nfinite, there is not a unique solution to the system. There are several 
methods of approach, but since the equations in each system are not all 
fequal weight, the A, and B, are assumed, as is necessary, to form a 
ecreasing sequence, and accordingly an iterative process is most appro- 

te. The method is described with reference to the system (32), but 
s applied in exactly the same manner to the system (27). 
The constants are evaluated on the basis of equality in the two series for 
values of 6. If the X. -series, which is convergent, is terminated at 


it is, at the value of m such that X,,,, = 0 for all practical pur- 


poses), then the A, -series must be terminated at ” s—1, and the con- 
stants evaluated by equating the two series, the final term in each being 
then coss@. The iterative method is to express A, (r > 0) in equations 


99 


32) in terms of A, and the known X,,, and then assign values to A, in such 


m? 
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a way that the two terminated series correspond exactly for values of § 
between 0 and 27. Thus, from equations (32), 
A, = 2/d,:[9(A9)—Ap9 Ap—(4+5) Xo] 
) F 2A, OF r 
A, 2/As 1} —Ap,, Ap—(a+b)X,- + {q(Ag)—Ap.p Ap—(a+b)X 


0,0 440 05 
se 


are substituted in the A, -series, A, being chosen to make the two series 
equal for @ = 4p, 6,,... and A, averaged over these values of 6. In practice 
the values of 6 are chosen to give points equally placed round the peri- 
meter of the cylinders as shown in Fig. 3. 


ay 











6, . 
Mg 
r=po(@) 
Fic. 3 
As an example, if the X,,-series is terminated at m — 3 (which is the 
minimum for any practical problem) and accordingly the A, -series at n = 2, 


the two expressions equated are 


c—(a+b)X,—(a- b)| X, cos 6 ; X, cos 264 X,, cos 36], 


and 
A o(Ago+Ao,1 COS 8) -+-2/A, [g(Ao)—Ap.9 Ag—(a- b) Xo |(Ay.9+-A;.1 COs A)cos 6— 
; 2A ‘ 
2 Az Nos A,—(a b)X, \ "9 fg(Ao) Xoo Ao (a+b)Xoj 
sis Fe 


(Ay 9+-As.1 Cos 8)cos 286. 


In practice, this method of evaluating the A,, is sufficiently accurate. The 
constants B, are similarly evaluated, the numerical work being slightly 
easier since the constant c does not appear in the equations. 

In view of the above analysis, it is easily seen that, although inclusion 
of further terms A,,,, (y > 1) in (1b) gives a more accurate representation 
of the Bessel function, the subsequent evaluation of the constants 
A,, B,, (n 0,1, 2,...) is correspondingly more difficult analytically. 

The stream function given by 

ob = ae”'?+-be-vP+ YY (A, cosnd+ B, sinné)K, (r/p), (33) 


n=0 
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vith the A, and B, (n 0, 1, 2,...) determined as accurately as desired by 


the above method, represents the flow past a general cylinder with line 


fsymmetry in the direction of the flow. The stream function satisfies 


the necessary condit ions 


i) V° |; /p* 
li) b> ae be iS 7 —> OO, 
ill) us c on the cylinder 7 poe(@), 
iv) Tj 0, where the circulation is taken around the perimeter of the 
ss-section of the cylinder. The velocity on the cylinder is given by q), 
here the g, and qg in equation (29) are obtained from equation (33). Stag- 
tion points on the cylinder are given by q 0. The stagnation stream- 
is given by , : 
‘ us ( g(a, b, X my ee A)), 
hich in the free stream near 2 ©, comes from ¥,, where 
Y.) p log lela {(i¢ a) h at |, (34) 


the given by equation (31). For numerical comparison with the results 
constant shear (5), we recall that in (5) the parameter V (= U/w/) is 


fined with / equal to y,,, where 2y,, is the maximum width of the cylinder 


rpendicular to the flow. The formulae in that case have y/y,, as dimension- 
ess length, so dimensionless quantities here must be multiplied by p/y,,. 


Ifa and 6 (the free stream constants) are given the values which result 


nthe free stream vorticity distribution being approximately a linear shear 


stribution, the stream function (33) gives an approximation, to any 
lesired degree of accuracy, to the requisite stream function for a general 


vlinder in such a flow. 
The case of the circular cylinder in such a flow, derived in a previous 
paper (1), is a particular case of the above method, with 4(@) As 


4. The general cylinder with line of symmetry perpendicular to the 

direction of the stream 

Flow with the type of free stream vorticity distribution given by equation 
5), past a general cylinder with its line of symmetry perpendicular to the 
lrection of the flow, can be obtained very quickly by taking the cylinder 
sthat given by (19) illustrated in Fig. 1. Since the y-axis need not coincide 
with the line of maximum width of the cylinder, a simple replacement of 
by x in equation (8) will not, in general, be sufficient. In Fig. 1 let the line 
imaximum width be at a distance 2, from the axis. The vorticity given by 


equation 8) is to be replaced by 


p | ae Lo+r)/P 1 he (o+z)P], (35) 
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where p is again taken as the most convenient typical length for the par. 
ticular cylinder under consideration, and a and 6b are the free stream 
constants. Thus, the vorticity in the free stream on the line of maximun 


width is —(a+5)/p? as before. Define the number s by s = 2,/p. Then 
from (35), a » 2 Act!» +. Be-* P] (36 
where A ae’, and B = be’. Fig. 4 illustrates the problem. 

| 





Y= + oo 
(Velocity profile) | { 





eS 








r=p@(0) 


The procedure for obtaining the stream function is that described in 
section 3, with A and B defined above taking the place of a and | 
respectively. In this case the Fourier expansion of exp{-+-¢4(@)cos 6} is 
needed in place of (24) and the stream function is of the form 

yy = Ae™?+ Be-*P+- JF A, cosnd K,(r/p), (37 
n 0 
where the A, (n 0,1,2,...) are obtained as before. All the physical 
features of the flow can be determined, and the general remarks made in 
the last section are applicable. 


5. The general cylinder in the ‘boundary layer’ shear distribution 
The general cylinder is again given by equation (19), and the function 
is introduced by the transformation 
% = acosy/l+bsiny/I+Y¥, 
where, from equation (13), ‘V’ satisfies 
i ! Al yl (38 
Cx" cy- 
and +0 as «+ —o. Fig. 5 illustrates the problem with a and b the 
values resulting in the approximation to the ‘boundary layer’ sheal 
distribution in the free stream. Proceeding as in section 3, a solution ot 
(38) is 
Y= Byé6Yi(r/l)+ > (A, cosné+B, sin nO)Y, (r/l), 


n=0 
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shere A,, and B, (n 0, 1, 2,...) are constants, and Y,, is the Bessel function 
fthe second kind. It is assumed that J, (r//) does not appear in the solution 
fom consideration of the case of the circular cylinder 4(6) 1, where 
eee - 
! uy 
a 
—_—__—_—_—» 
v 
r O + 
Yo r=pp(0) 








Xr=- 
(Velocity profile) 


lusion of J, (r//) leads to a contradiction of the independence of the Bessel 
inctions. The equation co1 responding to (23) is then 
a COS} p 14(@)sin GO) b sin} Pp 14(@)sin 6G} + B,@ Yip 16(6)} + 
> (A,, cos n6-+ B, sin n6)¥,{p/1d(9)}, (39) 
n=0 
here c is the constant value taken by % on the cylinder. The constants 
and B, (n = 0,1, 2,...) are obtained by equating the coefficients of 
nné and cosn@ as before. As in Appendix (a), the expansions of the sine 
nd cosine terms are obtained by expanding the functions in the form: 


ri 


) “ l J 1)?" . 9 
cos)” 4(8)sin 8 S Ram AS ) {b(A)sin 6}?", 
r=0 ual 
: " 1)" ])2r+1 : 
nd sin!? 4(6)sin 6! % See {d(8)sin 6}?r+1, 
\ 7 Pe (2r+-1)! 
T 0 


e Bessel function Y,{p//d(@)} is obtained as a cosine series using a com- 

able method to that in Appendix (b), resulting in a similar equation 
5b). 

Thus, following the method of section 3, the required stream function is 


btained, and is of the form 


acos y/l+-bsin y/I4 S (A,, cosnO+ B, sin n@)Y, (r/l), (40) 
0 


n 


and B, are known, and are such that »% satisfies 
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(ii) 6 > acosy/l+bsiny/l as r> o, 

(iii) & = c on the cylinder r = pd(@), 

(iv) T, = 0, where the circulation is taken round the perimeter of the 
cylinder. 

The velocity at any point in the flow can be obtained, and in particular on 
the cylinder by q,, defined by equation (29). Ifa and b are given the values 
obtained from (18), and if the length / used above is taken equal to (2/7) 
the field of flow represented by (40) gives an approximation to the flow past 
such a cylinder when the flow due to the presence of the cylinder is that ir 
an unretarded incompressible boundary layer of thickness 6, provided 
that, in the case when the axis of the cross-section of the cylinder is in the 
direction of the flow, (i) the length 2y,, is small compared with 64, and 
(ii) ¥)/d has a value not too near the end values 0 and 1 of its allowable range 

Davies (6) has recently sought to determine the magnitude of the dis. 
placement of the stagnation streamline in flows by traverses across boun- 
dary layers on a flat plate and a cone in sub- and supersonic flows. He used 
pitot tubes with circular cross-sections which he suggested must have 
radius less than 0-16. Thus, such a condition on p// gives 

Yb < 0-057. (41 

The corresponding deflexion y,, in the free stream, of the stagnatio 

streamline is given by equation (14) in the form 

c acos(y,/l)+-b sin(y,/1), (42 
where c is obtained from the circulation consideration, and a and } ar 
replaced by the appropriate expressions for the approximation to the shear 
flow in an incompressible boundary layer. Use of condition (41) in the 
numerical results is reasonable for comparison with Davies, since th 
defiexion in the free stream is in each case a result of the vorticity present 
in the flow. 

The flow past a cylinder with the line of symmetry of its cross-sectio! 
perpendicular to the direction of the flow, can be obtained in a similal 
manner to that described in section 4, with the appropriate changes in th 
flow in the free stream at infinity. For means of comparison, equation (4! 
must be altered so that the maximum length on the axis conforms with thi 
condition given by Davies. 

6. The elliptic cylinder r = p/(1-+-e« cos 6) 

The ellipse has p as the semi-latus rectum and « as the eccentricity. ‘Thu 

$(#) = (1+ecos@), from which it is easily seen that 0 < « <1 is 


necessary limitation on e. The practical difficulties which arise are in th! 


evaluation of the series (24) and (25). The appendix shows how this» 


carried out in this particular case. By using the appropriate equations, t! 
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fows can be obtained when the major axis of the cross-section of the 


ylinder is perpendicular to and in the direction of the flow respectively, 


where the flow may be that with the ‘linear’ shear distribution or the 


boundary layer’ distribution in the free stream. 
This particular case could be evaluated using elliptic coordinates and 
\athieu functions, but it serves asa practical example of the general method 


given above. 


7. The general cylinder of the form 7 = p(l+-y cos #) where |y| < I 

In this case, 6(@) = (1+ycos@). The appropriate equations (24) will 
be the product of the two series 

exp(sin@) and exp(ycos@sin@) {= exp(y/2)sin 26} 

expressed in their Bessel function expansion form, and (25) is the direct 
pplication of the method in the reference to Gray and Mathews (8) made in 
Appendix (6). 

This particular example illustrates the use of the general method when 
the cylinder has no natural coordinates. 

In each of the above examples implicit relations can be obtained for the 


deflexions of the stream lines and stagnation points on the cylinder. 


8. Discussion of results 

A method is derived whereby the stream functions for various shear 
flows can be obtained for any cylinder whose cross-sectional equation can 
be written in the form r pd(@), where 4(@) must satisfy certain conditions. 
The free stream flows are in terms of the general exponential and sine and 
cosine distributions of vorticity. By suitable choice of the constants in- 
volved, these can be favourably compared with a linear shear distribution 
nd that shear present in an unretarded boundary layer. 

The method is discussed in detail and the assumptions necessary for a 
solution to be uniquely defined are outlined. One important feature of 
ll the flows is that the stagnation stream line is displaced in the free stream 
irom the point on the axis opposite the line of symmetry of the cylinder. 

The two particular examples considered are those of the elliptic cylinder 
and the cylinder with equation of the form r= p(1+ycos@), where 

| is necessary to exclude singularities of the Bessel functions. In 
comparison with the former example, the latter is quite straightforward. 
The analytical difficulties in the elliptic cylinder are considerable, and the 
detailed analysis is given in the appendix. As shown in section 4, with only 
slight change in notation from the flow with the major axis in the direction 
of flow, the stream function can be obtained very quickly when the major 
axis is perpendicular to the flow direction. 
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If comparison is made with the experimental results obtained from 
pitot traverses in the boundary layers on flat plates and cones in supersonic 
and subsonic flows, certain restrictions on the dimensions of the cylinders 
in comparison with the ‘boundary layer thickness’ must be imposed. 

The examples cited do not exhaust the possible cross-sections for which 
the stream functions can be obtained using the general method described in 
this paper. They do, however, illustrate the difficulties involved in obtain- 
ing stream functions for non-uniform shear flows past cylinders whose 
cross-sections are other than circular. Further, it illustrates the difficulties 
which may be expected to arise in attempting to derive analytically the 
stream functions for variable shear flows of a more complicated nature, 


either compressible or incompressible. 
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APPENDIX 
(a) The Fourier ¢ rpansion of exp| sin@/(1 e cos 6)] 


The exponential function is 


; *(+1/).. 
exp|[ + sin 6/(1-+-e€ cos @)] ] S - [sin 0/(1-+-€cos@)]". (la 
rl 


Each term in the summation is expressed as a Fourier series, and these are then 
summed over 7 in the form given by (la). 
The expression (1-++¢cos@)-* (r 1, 2,...) can be expanded as a half-range cosin 
series in the form ‘ 
(1+-ecos@)* A, o/2 > a,,cosn, (2a 
n=1 


rn 


where, using Integral Tafel (7), a,,, is given by 


r—1 
; = WN [r+n—1\ (2r—y—2\ (1 —a*\” ‘ 
Ayn 2(1—a?)!—2"(1+-a?)"a2r+n-2 S | )\ : }( —} ’ (3a 
2 / y r l ] a* 
y=0 
where a lle (1 €? 1)?, (4a 


Ga . . ‘i 
and | } are the binomial coefficients. 
s 


The function [sin #/(1-+-€cos@)]" is an odd or even function of 6 according as r is an 
odd or an even integer respectively. If 7 is even, the function is expanded as a half: 
range cosine series, and if 7 is odd, as a half-range sine series. 

r even 
The function is an even function and can thus be expressed as 


[sin 8/(1+-€ cos @)]" sin’(a,9 2+ Sa 


~_ rn 


cos nO) 


A,,,cosm, {oa 
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| f ere = 
( ron r 
A. (2/zr) | [sin 8@/(1+-e cos 6)]’ cos mé dé 
eTSON] , 
0 
vlinders - 
: y 
dd. (2/ar) | sin’@{a, 9/2 = a,,,.co8 8) cos m6 dO, 
0 n=1 
r whicl 
+ . 7 xt 7 
x] iy P “ ys 
1bed | 7) | a, 9sin"@cosm@ d0+-(2/7) > a,, | sin"cosn9cosmé dé. 
obt ul 0 n=1 0 
whos fter a little manipulation, these integrals are evaluated using (7) to give 
iculties , 
; ‘ ( | - 
lly th , Sa > 
. , és 5 (7 m nv) 
nature 1=0 
™% 7 
a, { } m—n ” 
, a 2 }(r—m-+n) (if m n) 
1 
0 m—n>? 
matics 
7 7 yi 
ne pre ‘ 2 r ' 
S a ( ) ifn > m) (6a) 
oud ys i(r—n-+-m) 
7 1 
it must be remembered that r, m-+-n, m—n are even integers, and a,,, for all 
01s gi by (3a). In any calculation, r and m are given, and the varying ” in 
DON I is chosen so that m nm, mm mare even. 
function in this cas un odd function and is expressed as 
ire tl sing €cos@)|" ae Bes sin mé. (7a) 
m=1 
cosir roceeding as above, the coeffi nts B are given by 
Tr 7 
2 ; 2 7 
B S a | 
— Z 3( m rn) 
i i) 
l ae 
2 7 
) ‘*\ a ( 1 n 7” 
, Fs (r m n) (if m n) 
? 1 
f 0 m—n>? 
ti 7 


+ a r } ifm < n) (Sa) 


r, m+n, m—n, are odd integers. 


Using equations (5a) and (7a), the required Fourier expansion of (la) is now 


i n as 


exp n@/(l ecos? 


Law (2r)! Lat Loe (27 +1) 
m=1 r=0 m=1 


aie a “ am Ao, oie <A B., ' 
.. oo > =m cos mé + S * hm sin mé, (9a) 














422 J. D. MURRAY 


where A,,, and B,,,, for allr and m are given by equations (6a) and (8 a) respectively, 
It should be noted that (9a) is only valid for 0 e < l. For convenience, defining 
X,, and Y,, by the following expression, equation (9a) is written in the form 
x x 
exp[-tsin 6/(1+ cos 6@)] Y X,cosmO+ > Y,,sinmé. (10a 
m=0 m=1 
These summations converge strongly, with stronger convergence for smaller values 
of «. The numerical values of the coefficients X,, and Y,, are easily calculated: as an 
example, when € 0-25, 
exp| L_sin 6/(1 0-25 cos 6) j 1-2799 —0-0727 cos @— 00-2715 cos 204 
0-0732 cos 38— 0-0088 cos 46...+- 1-1495 sin 6+ 0-1630 sin 20 + 0-0210 sin 36... 
(b) The Bessel function K,{1/(1+-€ecos@)] as a cosine series. 


To a first approximat ion K,,( 1) may be used, since the argument has limits 1/(1 Le), 
Thus, for small e (e < 0-1) K, may be considered constant, and the method outlined in 
(1) is directly applicable. The discrepancy from the true value increases with € and n, 
Writing. K,[1/(1+¢cos@)] in the form 
x 
K,[1 (] ecos@)} em An, cosy), (lb) 
y=0 

the procedure for determining A,,,, as a function of € is carried out as follows, using 
Gray and Mathews’ Treatise on Bessel Funciions (8). The expansion of K,(R) asa 


cosine series is given as 


K,(R) (2R/(ab))"(n—1)! ¥ (n+ p)C#(cosO)K,, (b)Iy. (a), (2b 
p—0 
where R? a®?+-b?— 2abecos6@, I is the Bessel function of the first kind, n + 0, 
a b|, and C!\(cos @) is Gegenbauer’s polynomial, which for integral p and n may be 
defined by 
1 (Dos)P 9 ) « 
CX) (n+p 1)! (2v) [1 _P(P 1) (2v)-2 P(p- 1)(p <)(7 3) (2v)-4 
: (n—1)! p! (n+p—1) 2!(n+-p—1)(n+p—2) 
from which C(v) l, and C7(v) 2nv. The case n 0 is given by 
K,(R) Ii(a)Ky(b)+2 S K,(6)1,(a)cos pd, (3b 
p=1 
where |a b|. Comparing these formulae with the Bessel function in question 
ne (1+ecos@)-? kay 9+d,,c0s 8, 


to a first approximation from (2a), with a, and a, , given by (3a). This gives reason- 


ably good results if € 0-5, the smaller the value of € the better is the approximation 
Using the last relation, comparison with R? in (2b) and (3b) gives 

a? +b? 4-9 2ab A515 
which gives a 4[(442,9— 42,1)! —(402,9+4,,)#] 
and b = 3(44o9—a2)!+ ($a2,9+2,)!), (4b 
since |a| < |b|. If a and 6b are to be real, it ean be easily shown from (3a) and (44 


that e€ must be less than 0-5. For larger values than 0-5 the first approximations givé 
poor results and the values must be suitably amended by trial and error. In equation 
2b) R appears in the form R", and so further use of (2a) must be made, the cast 
ss 1 being R 1/(1+€cos@) ka, 9+a,,cos 8. 

Substituting the required quantities in (2b) and (3b), approximate expressions for 
A, are obtained. Table 1 (6) illustrates the form of A,,, and A,, , for various values of n, 
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ctiveh where a and b are given by (4b). Terms with y > 1 are similarly obtained from 
sation (2b). 
TABLE 1 (bd) 
0 ' Ano 
o\9) 
I 1b)[ ay 9 1,(a)K,(b)+ 2a, , 1,(a)K4(5)] 
lp 9 1(a)K(b) + 6a, (a) K,(0)] 
Ana 
‘ 1\9) 
. I ay 1 1,(a)Ky(b) + 444 9 12(a)KQ(d)] 
b)*[ 12g 9 13(a)K3(b) + 2a91 J(a)K4(b)] 
n tl uluation of the constants A, and B,, in equation (23), certain analytical 
ilt rise as a result of K,[1/(1+-ecos8@)] being itself a function of 6. Accord- 
ngly, if the series (1b) is terminated at y 1, A,» and A,,, must be amended by 
| and error to give a reasonable approximation to K,[1/(1-+-e€cos@)]. The pro- 
lu to evaluate A,,, and A,,, from the formulae given above, and to amend the 
rica ilues for each » and e. Table 2(b) illustrates the numerical changes 
R necessarv in the case of K 1/(] 0-25 cos @)]. 
TABLE 2 (b) 
«| { An,1 
‘ e? T heore lmended 
é lue Mes 
4227 I 3 0°1603 
5978 o-I174 0°2424 
I°7IO! 0*7004 | 03954 
Figs. 1(b), 2 (6) illustrate the graphs of K,[1/(1+ecos@)] with e = 0-25 and n 0, 1 
9 respectl\ Vv, with the various approximations to the exact curve. 
6 - 
' | 
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A more accurate representation of K,[1/(1+-€cos6)] is obtained if A, (y > 1) are 
included in the summation, but with a proportional increase in analytical labour, 










































It is sufficient, however, for the purposes of illustrating the method and analytical 
08 T 
Amended 
a approximation 
*s) L Exact 
™ 
‘Y 
“—- 0 7r \ 
r ~~. 
|O ss 
o “a 
-_ Q inerncrtartmer) ‘1. 
O approximation™, 
“ 
035 47 87 7 107 147 27 
9 9 9 9 
A 
Fic. 2 (b) 


difficulties involved to assume that the series given by (1b) terminates at y = 1, 
y j Oo r ~ 

etving k,[1/(1+€ cos 0) | Ano An cos 8, (5b) 
where A,,) and A,,, are the amended values giving the best approximation to the Bessel 
function for a given e. Equation (5b) gives a good approximation without amend- 
ment to Phi for small values of e, the accuracy decreasing with increasing e. The par- 
ticular case of « = 0 (that of the circular cylinder), gives, as may be easily checked 


A , 0 for y > 1. 
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THE STEADY TWO-DIMENSIONAL FLOW OF 
VISCOUS FLUID AT LOW REYNOLDS NUMBERS 
PASSING THROUGH AN INFINITE ROW OF EQUAL 
PARALLEL CIRCULAR CYLINDERS 


jy K. TAMADA 
(Department of Aeronautics, University of Kyoto, Kyoto, Japan) 
and H. FUJIKAWA 
(Jr. C. of Eng., University of Osaka Pref., Osaka, Japan) 
Received 10 May 1956; revise received 21 February 1957] 


SUMMARY 


\ detailed discussion based upon Oseen’s equations of motion has been made of 
steady two-dimensional flow of a viscous fluid passing perpendicularly through 
ninfinite row of equal parallel circular cylinders regularly spaced. It has thus been 
nd that the drag acting on any one of the cylinders in the row is always greater 
n that acting on the same cylinder when it is immersed alone in an unlimited 
niform flow with the same velocity. In particular, when the Reynolds number of 
flow is sufficiently small, the drag is found to be proportional to the flow speed U, 
le it is proportional to U/log U in the case of a single cylinder. Thus, the inter- 
rence effect between the cylinders in the row is very remarkable at low Reynolds 


1, Introduction 

Tue steady flow of a viscous fluid past a cylindrical obstacle at small 
Reynolds numbers has hitherto been studied theoretically by many authors 
on the basis of Oseen’s linearized equations of motion. However, so far 
is the present writers are aware, most of these investigations concern the 
case of flow around a single obstacle. 

It is natural to expect that the disturbance due to an obstacle will spread 
very far in the flow field at low Reynolds numbers (e.g. the flow of very 
viscous fluid). Therefore, when two or more obstacles are present in the 
otherwise uniform stream, the interference effect between them will be 
very conspicuous at low Reynolds numbers. In this respect, the present 
writers have studied the steady two-dimensional flow of incompressible 
viscous fluid past an infinite row of equal regularly spaced circular cylinders 
ying in a plane perpendicular to the uniform stream. The flow of this type 
seems to be of considerable interest both from theoretical and practical 
points of view. The whole analysis is based on the Oseen approximation. 
Denoting by R and R,, the Reynolds numbers referred to the diameter of 
each cylinder and the central distance between the two neighbouring 
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cylinders respectively, we may distinguish the case in which R is small 
but R, is not small from the case where R and R, are both small. In the 
former case, a solution in the form of R-expansion is obtained correct to the 
order of R?, the lowest order term being O(1). In the latter case, the solution 
is expressed in the form of power series in the diameter-distance ratio of 
the row. 

Detailed numerical calculations are made for the drag coefficient at 
various Reynolds numbers ranging from 0 to 8, and for values of }, 5, 3. 


=\5, for the diameter-distance ratio of the row. 


2. Fundamental equations and boundary conditions 

Let us consider the steady two-dimensional flow of a viscous fluid passing 
through an infinite row of equal circular cylinders of radius a, spaced 
regularly at a distance h between the centres and lying in a plane per- 
pendicular to the uniform stream. We take the rectangular coordinates 
(x,y) in the plane of fluid motion, in such a way that the origin coincides 
with the centre of any cylinder of the row, and the x-axis is along the direc- 
tion of the uniform stream. 

Let U be the undisturbed velocity, u, v the x, y components of the 
velocity at any point, and ¢ the vorticity defined as 


y 


C (év/ex)—(eu/ey) 27(ew/0eZ) F 
w u—twW, 2 uti, Z=x—wl " 
Then Oseen’s equation of motion may be written as 
(V2—2ke/ex)l = 0, k = U/2v, (2 


where v is the kinematic viscosity. 
From these equations it will be seen that the complex velocity w for the 
present problem may be put in the form: 


where 
W,= > A,,(a/z,)” 


m=1 
WD ekr,cos 8, ~~ a_K (ker Je imé, (? (3a 
s pa mAm\FTs : 
U0 x. 
~ reo. Ze z—ish = r. ids 


and ¢, A,,, and a, 


m? v 


are arbitrary real constants to be determined by the 


boundary conditions and K,,,(kr,) denotes the modified Bessel function. 


Further, since the vorticity ¢ should be real, re(@w/éZ) = 0, whence 


m 


9 
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Now, for |x|» 0, w should be equal to the undisturbed velocity U. 
from (3) it will be seen that when 2 > —« 


w~ U+e > (A, a/z,) ~ U+e—(na/h)A,. 


Hence € (7a h)A,. (5) 
\lso, it can be shown after some calculation that 


~ 


w—> l (27a /h)A,+-(27/kh) > An 
m==0 
for «> +-00. Therefore the following relation is expected to hold: 
A,+(ka)- > A», 0. (6) 


The relations (3) together with (5) and (6) already satisfy the condition 
t infinity as well as the condition of periodicity in the y-direction. The 
nly remaining condition to be satisfied is that w = 0 on a cylinder r = a. 
For the fulfilment of this condition, we expand w in a Fourier series} in @ 
nr =a. Then we equate every coefficient of exp(in®@) (n 0, +1, +2.,...) 


zero. Thus, after some calculations, we get the following simultaneous 





juations 
L+t'°s g Giy—n))4m = 9 (n = 1,2,...) (7) 
l (n 0) 
VY FA S G (3) 
1 ee ‘Jmn mn Mm 10 (n oe i 
where F 0 7a h 
l : B (=) : ( ven) 
( ‘Sides 1 | - s8--n = eve 
F | sini(st+n) *™ ss h 
| 0 (stn = odd > 3) 
q K,, (ka)L,, (ka) (9) 
G 2 SF (—1) GMT, (kh) I (ka)h,,,,(ka) 
| 0 (l+-m = odd) 


ind B, denotes Bernoulli’s number. 


t For this purpose, it is convenient to use the addition theorem for K,,(kr,)exp(in9,). See, 


for example, G. N. Watson, J'heory of Bessel Functions, p. 361. 
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The first equation of (7) verifies (6), as expected. Elimination of the 4 
from (7) and (8) gives 


2 —U (n= 0), 
> Xmn An ( hy (10 
m= | 0 (s = I, 2....), 
x 
r re y, a y — ; i" s + y t € 
W here Xmn Imn ' | Zz FAI m(-s) | Gin(-s)f- (10a 
s=1 


Also, the a,, must satisfy the equation (4). Thus, equations (4), (7), and 
(10) determine in principle the constants a,,, A,, and hencethe solution of the 
present problem. It is, however, difficult in practice to solve these equa- 
tions exactly. On the other hand, Oseen’s equation itself is valid only at 
small Reynolds numbers, and our main interest in the present work lies 
also in such flows. In what follows we accordingly seek for approximate 
solutions valid for the flow at low Reynolds number. 


3. Approximate solutions for small Reynolds numbers 

There are two parameters ka and kh in the present problem, and we 
shall first treat the case in which ka is very small but kh is not. In this case 
we expand the coefficients y,,,,, in (10) in powers of ka’ by expanding J, (ka) 


and K,,(ka). Then it is found that a,, = O[(ka)?”"| and the unknowns 


a,, can be determined successively in the form of power series in ka. Thus, 


after simple but tedious calculations, we obtain the following result: 


(a,/U) (2/Z)+-(2/Z*)(Q?— QZ-+-(§)(2a/kh)?+- 1} (ka/2)?4 of 
(a,/U) = (2/Z)(Z—Q)(ka/2)?4 | (ll 
where Z = 4T,—(2n/kh)+28+1 
Q 477,-4-(1)(2r/kth)®—(2ar/kh) +(4) 
(lla) 


Tas 2, Kaslskh) (s= 6 8...) 
S = log(2/ka)—y (y: Euler’s constant) 


Further, (7) combined with (11) determine the constants A,,. But we shall 
only note here that A O| (ka)m-*}. 
We now proceed to the evaluation of the drag D per unit length of a 


m 


cylinder in the row. By summing up the stress acting on the surface of a 
cylinder, it can be shown without difficulty that 


D = 2npaU A,. (12) 
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Substituting A, from (6) and (11) we arrive at the result: 


; D Sa if > , o2[ am \? \ Rk? é 
( nD - ] Z Q)*4 : ie OS = pee 
2apU2 RZ Z\ é) 3 (F) | 64 | | (13) 


R = 4ka = 2aU |v, R, = 2kh = hU/v 

ere Z and Q are givenin (lla). If we here take the limit as h > «0, we 
btain a formula for the case of a single cylinder in an unbounded stream, 
shich is seen to agree, up to the order considered, with the result given 
y Tomotika and Aoi (1). 

Further, it can be shown that the pressure drop in the far down-stream 
from the row of cylinders is given by the formula 

(p Pp ~)/pU* (2a h)Cp, 

where p_.. and p ‘e the pressures at x 1 and x |-oo respectively. 

We next treat the case in which ka and kh are both small. We begin as 
in (10). 


First, it is necessary to get asymptotic expansions for the function 


before with an approximate evaluation of the coefficients y 


I,(kh) for small values of kh. To this end we rewrite T,,(A) in the form 


x 


’ = i: £ 
T (A) T'_,,(A) > K, (8A) = | 15 7) | 


l —_ 





1 
exp l(t -) | #2i-1 dt, (14) 


Expanding the integrand in powers of A, carrying out term-wise integration, 
ind summing up the resulting series, we get the expressions 
Ty(A) ~» (4)(427/A) + (4)log(A/427) + (y/2) 
—£(3)(A/42r)2+- O(A!) - 
T,(A) ~ (2,)(42r/A)2—(4)(427/A) +. (4) + O12) , (15) 
T;(A) ~ (B;/8j)(427/A)?! —(B;_,/4)(427/A)?/-2+- O(A-* +4) 
(j = @, 3...) 


where C(7) is the zeta function and the B, are Bernoulli’s numbers as before. 





3y making use of these results together with the expansions of J,,(ka) 


m 


and K,,,(ka) we can obtain, after tedious calculations, required expressions 


ior the y,,,,,.. However, we shall only note here the following results: 
X2m,21 O| (ka ) ass |, X2m+4 1,2” Of (ka) 
X2m,2n 1 O| (ka ) sash |, X2m +1,2n+1 O| (ka) — 1 ws (16) 
(m ie Dinan © (n O, 15052) 


From this and (10) it is suggested that 
As O| (ka)?™], Panas O[ (ka)? +2}, 


m 


therefore, to the lowest order in ka, we can determine the a,,, separately 
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from the ay,,,,. Thus, after some calculation, we obtain 
ay/U = —(2/Ay)+O[(ka)?*], dy/U = (Ag/Ay)(ka)?+- Ol (ka)*] ) 
a,/U —(A4/Ag)(ka)*+ Of (ka)®], ... J 


where 


» (li 





Ay 1—2logr+i7?—,h7*+- ads 7® 
3156007” I 54432007 1 O(7"*) 
Ay is(t” ‘ I ‘i reas) +t O(7'°) be (17a 
Ng se0n0(7* 217°) + O(7*), 
T 2r7a/h 
Inserting these in (12) combined with (6), we arrive at a formula 
D/(pU) RC) 87/A,+O0(R?), (18 


where yu is the viscosity coefficient, A, is given by (17a), and 


R tha 2aU /v 
as before. 


4. Numerical discussion and conclusions 

We are now able to calculate the drag per unit length of a cylinder in 
various rows with different diameter-distance ratios 2a/h, at arbitrary 
Reynolds numbers &. For sufficiently small R, formula (18) is used and 
the formula (13) is valid when R becomes comparable with or exceeds the 
ratio 2a/h. It may be remarked here that (13) diverges for some value of 
R which makes the denominator Z vanish. This is, however, a trivial 
singularity brought about by the expansion of the modified Bessel function 
K,, (ka), and is of common occurrence in formulae of the same kind (ef. 
reference 1). To avoid this difficulty it is better to use a direct method 
in which we start with the simultaneous equations (10), evaluate their 
coefficients and then solve the equations with numerical coefficients thus 
obtained. 

Tables 1-4 and Fig. 1 show the variation of D/uU’ with the Reynolds 
number RF for various values of 2a/h, while Table 5 gives the values of D/uU 
against 2a/h for sufficiently small R for which D/uU becomes independent 
of R. In thesetables (D/uU) 9, (D/uU),, and (D/nU), denote respectively the 
values calculated by (18), (13), and the direct method mentioned above. 

These results may be summarized as follows: 

(1) The drag acting on a cylinder in the row is, in general, greater than 
that on the same cylinder when alone in a uniform flow of the same velocity. 
Further, the drag increases as the diameter-distance ratio of the row 
becomes large. 


(2) As the Reynolds number of the flow is reduced, the value of D/uU 
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deviates more and more from that for the case of a single cylinder. More. 
over, it tends to a finite limit when R -> 0, whereas in the case of a single 
cylinder it tends to zero as (log R)-}, 

Lastly, the following remark may be made. The flow is narrowed down 
as it passes through the row and the shear in the fluid near the surface of 
each cylinder thereby becomes greater. This brings about the increase in 
drag on each cylinder. In this respect, there is an essential difference 
between the case when the number of cylinders is infinite and that i: 
which it is large but finite. In the latter case, the flow may pass round the 
system of cylinders, leaving the fluid inside the system almost at rest. Thus, 
in such a case the drag on a cylinder inside the system is expected to be 
very small. 

The writers wish to express their cordial thanks to Professor S. Tomotika 
for his continual interest and encouragement throughout this work. 
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& singl ON THE EQUILIBRIUM OF A STRATIFIED 
by » 4 4 
LAYER OF FLUID 
d dow) 
ey: By B. R. MORTON 
rease (Department of Mathematics, University College, London)+ 
fferenc Received 16 October 1956] 
that i SUMMARY 
und the his note demonstrates the equal importance of viscosity and thermal conduc- 
t. Thus ty (or any equivalent process of diffusion) in determining the departure from 
1d to he librium of continuously stratified layers of fluid. The analysis is carried out by 
xtens n ol pre viou Val itional principles, and IS applied to freely bounded 
| layers which are thermally stratified ; the initial rate of growth of small 
is - 
Motika rbances is shown to depend on the Rayleigh and Prandtl numbers. The method 
‘kk. isual critical Rayleigh number for the onset of convection in unstably 
d layers with a linear density gradient, but it also indicates that this critical 
the Rayleigh number suffers little reduction when the stratification is 
ginally non-uniform. Although the analysis is applied to a system which cannot 





| physically, the corresponding results for real systems should be qualita- 


simular. 


1. Introduction 

. horizontal layer of fluid is heated from below it becomes unstable to 

ill disturbances when the temperature gradient is sufficiently large. 

e equilibrium of such thermally stratified layers of fluid was first studied 

Rayleigh (1); since then notable contributions have been made by 
effreys (2, 3), and by Pellew and Southwell (4) who have also given a 
survey of the earlier work. In the theoretical treatment which has been 
leveloped it is assumed that there is a linear temperature profile in the 

ver, and critical conditions are determined under which a small dis- 
turbance will first be amplified. It is found that the condition for marginal 
stability can be described by a constant value of the Rayleigh number, the 
particular value depending on the nature of the surfaces bounding the 

ver, 

\lthough the Rayleigh type of analysis has been used successfully to 
pedict the critical parameter for instability of the steady state, there are 
problems to which it cannot be extended readily. To overcome this diffi- 

ity Chandrasekhar (5) has derived a variational principle and has sug- 
gested that this could be used as the basis for an approximate treatment 
{the equilibrium of stratified viscous fluids. Subsequently Hide (6) has 
pplied this method to two problems in which layers of immiscible fluid 
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are superposed, and to a third case in which a layer of viscous fluid js 
stratified continuously. However, Chandrasekhar developed his variational 
formula for immiscible viscous fluids, and consequently it takes no account 
of any diffusion of the property causing variations in density. It cannot, 
therefore, be used in considering continuously stratified fluids, whether the 
stratification is produced thermally, by dissolved solids, or in any other 
way. 

The purpose of this note is to bring out the equal importance of the 
diffusive effects set up by gradients of velocity and of density in determining 
the departure from equilibrium of an unstably stratified layer. To illustrate 
the general nature of the equilibrium it will be sufficient to consider the 
stability of a freely bounded horizontal layer under various conditions of 
thermal stratification and with slightly simplified boundary conditions, as 
the effect should depend primarily on the physical properties of real fluids, 
Moreover, it will be shown that for this limited application it is relatively 
easy to extend the previous variational treatments, given by Rayleigh (7 
for inviscid fluids and Chandrasekhar (5) for viscous fluids, to cover the case 
of real fluids possessing properties of viscosity and thermal conductivity, 
The results of this analysis will show that the initial growth of small dis- 
turbances in continuously stratified layers is characterized by both the 
Rayleigh number and the Prandtl number (i.e. the ratio of kinematic 
viscosity to thermometric conductivity), and not by the Grashof number 
alone as Hide has predicted. It will be shown also that a non-linear ten- 
perature profile can have little relative effect in modifying the critical 
value of the Rayleigh number for marginal stability. 


2. The equations representing the motion 

Consider the equilibrium of an infinite horizontal region of incompressibl 
fluid of depth h, with its upper and lower surfaces free. The system will bé 
referred to rectangular cartesian coordinates (a;) = (x, y,z) with the origin 
in the mid-plane of the layer and the z-axis directed vertically upwards 
and it will be assumed that the acceleration due to gravity (g) is constant 
throughout the layer and is directed vertically downwards. It will be 
assumed also that the coefficient of viscosity » and the thermal conductivity 
are constants, and that local variations in temperature are always small it 
relation to the absolute temperature of the upper surface (and hence that 
variations in density are small relative to p,, the density near the upper 
surface). 

The equations representing the convective motion set up by a disturbance 
in the fluid can now be written in the approximate form (see, for example 
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Pellew and Southwell (4)). 





cu ou C 2. = 
Pizt ii; APo™ P). A; pg+pVu; 
ct CX Cx: 
J D) 
cu 
0 
On : (1) 
o(T +-@) o(T' +0 ora 
La ) — eV T+) 
ct CA 
p = po(1—Pf8) 


where V? 


the undisturbed, and 7’ 


T(z) 


# the disturbed, absolute temperature, py the 


C7 /ex,Ox,;, (u (u,v, w) is the disturbance velocity, 7 
undisturbed and p the disturbance pressure, p, the undisturbed and p the 
listurbed density, A = (A,) the unit vector along the z-axis, 8 the coefficient 
f cubical expansion, « the thermometric conductivity, and t the time 
variable. 

The previous analyses have been based on equations (1) and have been 
oncerned with the equilibrium of a stratified layer of fluid which is initially 
ina steady state, and thus in each case the undisturbed temperature profile 
has been taken as linear. Up to the present the analysis has not been 
extended to determine the behaviour of layers in which the stratification 
of density is continuous but not uniform. According to the above assump- 
tions, non-linear temperature profiles cannot represent equilibrium con- 
litions and, therefore, must change steadily with time. However, it may 
still be possible to find the effect of small disturbances on these non-uniform 
iyers, provided that the time scale for the decay of the temperature profile 
isassumed to be large in comparison with that for the initial growth of the 
listurbances. Some of the consequences of this condition are considered in 
the appendix, but it may be remarked here that this type of analysis applies 
nly to infinitesimal disturbances, and ceases to be valid as soon as there is 
significant growth of the motion. 


In the undisturbed state (w;) and @ are zero, so that equations (1) reduce to 


Po 


CL; 
i 


A; Pog 
wey 


P Po 
In the disturbed state assume that (w;) and 6 are small quantities of the 


rst order: their squares and products are then of the second order of small- 


hess and may be neglected in relation to first order terms. From equations 
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(1) and (2) it follows that 


7 
(: = VV? \u, — <P + Boa; 0, (3 
ot Py OX; 
et (4) 
Ca P 
¢ \ ial 
(' _ v9 = —w@t. (5) 
ot dz 


where v is the kinematic viscosity, and the undisturbed temperature T js 
taken to vary only with height. Taking the divergence of equation (3), it 
follows from (4) that 
Es fale] 
0 V?p-+ Bg; - 


Py Cx 


| 


and hence, operating with V? on equation (3), 


l¢ ‘ 2A \ 
pV) Vu; Bglr, V20—A, —— —). (6) 
( ; aii, ” 6x; Ox; 
Then, from (5) and (6) 
(: KV? (- : vV2) Vw » —fig =” Vite (7 
ct ct ~ dz 
where V2 a i 


ew 9° 
Cx" Cc y 


In order to determine the initial rate of growth of a disturbance of given 
wave-number, equation (7) is reduced first to an ordinary differential equa- 
tion in z by separation of the vertical and horizontal variables. This can be 
done by assuming that 2 9 

ep) i Viw+ kw 0, (5 
and that w(x, y, z) is of the form f(x, y) w(z). Then in any pattern of motion 
which is such that each element of fluid moves within a restricted region, 
k must be real. For, from Green’s theorem, 

k? | | w* dady = — | | wV? w dady 

: a 
re LT few\? ow\? r Ow 
| | ( A: dudy > w— ds, 
Je OX cy ; y On, 
I II 
where I is the whole area of the plane considered and II the perimeter 0! 
this area. If the boundary is chosen to coincide with that of a set of con- 
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vection cells (of any kind) the final integral will vanish, as at each point 


nm the perimeter either the vertical velocity (w) or its normal gradient 

w/én,) will be zero: hence k? is positive, and k real (cf. Pellew and Southwell, 
4)). Thusk may be taken as the wave-number of the (harmonic) disturbance 
in the horizontal plane. This method of analysis is equivalent to reducing 
, general disturbance into a set of normal modes, each characterized by its 
wave-number k, and this may be done because the equations have been 
inearized. 

If the amplitude of the disturbance of wave-number k varies as exp(nt), 


equation (7) reduces to 
[ —n+«(D?—k?) || —n+-»(D?—k?*)|(D?—k?)w = yk®wDT, (9) 


where D = d/dz and ; Bg. If an exact solution is krown for w/z), 
equation (9) gives the corresponding values of n for each k. However, the 
exact form for the profile of vertical velocity is known only in the case of 
i linear temperature gradient, and so to obtain an approximation for » in 
the more general case considered here it is convenient to develop equation 
9) into a variational form. 

The boundary conditions which will be satisfied at the free bounding 
surfaces are 

(i) w = 0; this applies strictly for a rigid surface at which there is slip, 
but although the consequent small increase in constraint on the growing 
modes may be expected to modify the quantitative results of the analysis, 
the qualitative results should be unaffected ; 

ii) Cu/ez, Cv/ez = O (zero viscous stress), whence from (4) é?w/éz? = 0; 

(ili) @ 0, hence from equations (8) V?@ = 0, and from (6) 


(e/et—vV?)V2w = 0. 
Thus at a free boundary, 


w(z 0. D2w(z 0, D4w(z) 0, D&w(z) = 0, 


3, A variational principle for the layer with free surfaces 


Equation (9) can be written 


i 


n,+«(D*—k?*) || —n,;+v(D?—k?) |(D?—k?)w; yk?w, DT, (10) 


where the velocity w; corresponds with the characteristic value n,; of the 
growth factor in the exponential term. Multiply (10) by w; (corresponding 
with n,) and integrate over the range —4h < z < 4h; then, after repeated 


integration by parts, it follows that 


I, = 0, (11) 


tj 

















438 





B. R. MORTON 


where 


Lh 
I;; n? (Dw, Dw,+k*w, w;' dz 
“th 
lh 
(k-+v)n; | LD?w, D?w;+ 2k? Dw; Dw;+-k*w; w;} dz + 
Yh 
bh 
tev | {D§w, Dw;+3k?D?w, D?w,+ 
‘th 
+ 3k*Dw,; Dw; + k8w,; w+ (yk?/Kv) DT w, w;\ dz. 
Hence, subtracting J;; from J,,, if n; 4 n; it follows that 
bh 
(n;+n;) | {Dw, Dw,+kw,w;} dz 
‘yh 
bh 
(v-+-«) | (D2w, D?w,+ 2k? Dw; Dw;+- kw; w;} dz = 0. (12) 
a 


A complex value of » will be associated with a complex value of w; their 
respective conjugates (7, #) will also be associated, and will satisfy the same 
equations and boundary conditions. Thus if ”,;, w, are taken as n, w, and 


n;, Ww; as H, WH, respectively, equation (12) reduces to 


bh 
2 re(n) {| Dw)|?+-k?\w!? dz 
—th 
bh 
(vt+K) | {|D?w|?+ 2k?|Dw|?+-k4\w|?} dz = 0. (13) 
“th 


Hence if » is complex it follows that the real part of n must be negative; 
that is, instability through oscillations cannot occur because of damping. 
Therefore re(n) > Ois possible for real values of n only, and imaginary values 
of n need not be considered in a search for unstable modes. It also follows 
from (11) that if x is to have real positive values (and these give rise to the 
only instabilities) D7’ must be negative over at least part of the range. 


To obtain an equation which can be used as the basis for a variational 


determination of n put i = j in (11) and suppress subscripts, 
bh kh 
n? ( ((Dw)?+-k?w?} dz +-(«+v)n | {(D?w)?4+ 2k?(Dw)?+ k*w*} dz 
th th 
th 
+Kv | {(D§w)?+-3k?(D2w)?4-3k*( Dw)? + k’w?+- (yk?/«v) DT w*} dz = 9. 
—in (14) 


Equation (14) is in suitable form, for if it is written 


n*I,+-(x-+v)nI,+x«vl, = 0, 
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then to the first order, if dw is an arbitrary functional variation in w (com- 


patible with the boundary conditions on w), 
(2nI,+(«-+v)1,}dn+-{n?6F,+(k+v)n 61,+-«vdl,! = 0. (16) 


On integrating by parts in each case it is found that 


h 


SI, 2 | 8w(D2—k2)w dz, 
h 
h 
I, 2 | dw(D?—k?*)*w dz, 
bh 
h bh 
él, —2 du(D?—k* Bw dz +-2(yk?/«kv) | dwDTw dz, 
h bh 
nd hence 
I, + (k+v)nbh,+-«v dl. 


—? {«v( D2— k?)8w— («+ v)n( D?— k*)2w-+- n?(D? -k?)\w—yk? DT w} 


x bw dz. 


\s any function w which isa solution of the problem satisfies the differential 


uation 
xv( D?— k?)8w— (« +-v)n( D?—k?)?w+- n?(D?—k*)w—yk?DTw = 0, 


t follows that if w satisfies the differential equations representing the 


otion together with the boundary conditions for free surfaces, then 
(2nI,+(x«+v)L,} dn 0. (17) 


But the factor {2nJ,+(«-+v)J,} in (17) is a positive definite function under 
these conditions; hence 6n = 0 and vn is stationary. 

Ifthe correct solution for w is inserted in equation (14) the correct value 
iorn will be found. However, if the precise form of w is unknown it follows 
tom the stationary property that any reasonably close approximation will 
sive a good estimate for n from equation (14), since all small variations in w 
vhich are compatible with the boundary conditions leave n unchanged to 
he first order of accuracy. 

It is apparent from (14) that the effects of thermal conductivity and 

scosity enter symmetrically through the coefficients « and v, and that in 
general neither may be neglected in considering the equilibrium of con- 
tinuously stratified fluids. Although no allowance has been made for 
riations in « or v (as in Chandrasekhar’s formulation) this neglect should 
aly introduce small errors into the numerical results which will be obtained. 
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4. The approximate treatment of equilibrium 

3efore considering the various ways in which the equilibrium of the layer 
may be affected, it will be convenient to change the variational equation 
into dimensionless form: this can be done by means of the transformatior 


Z hZ, k . kK, n (xv)! N , wz)= WZ), D d ld ; 
h h? dz hdZ 

Two non-dimensional parameters can be used to describe the equilibrium 
and the growth of disturbances in the layer: the Rayleigh number 
A = yh’ AT'/xv provides an ‘entire’ parameter, and the ‘profile number 
b(z/h) = hDT/AT provides a local parameter, where AT’ T,—T,| and 
1), T, are the temperatures of the lower and upper surfaces respectively. 
A does not depend on local values in the temperature profile, while } 
represents the ratio of the local to the average temperature gradient. 
Equation (14) can now be reduced to the dimensionless form 
N? | {((DW)?+ K?W?} dZ-4 

+-(ot+o-*)N ( ((D?W)?+-2K?2(DW)?+ KAW?) dZ 


1 


| ((D3W)?+-3K2(D?W)?-+-3K4(DW)?+ K8W2+ K2AbW?) dZ = 0, (18 


where o = v/« is:‘the Prandtl number. 

In previous treatments no attempt has been made to determine the form 
of motion that may develop in unstably stratified layers where there is a 
non-uniform distribution of temperature. However, for the case of a linear 
profile Pellew and Southwell (4) have given exact solutions for free layers 
in which there is neutral stability (i.e. N = 0). They distinguish between 
‘even’ solutions W = C cos(2p+1)7Z, and ‘odd’ solutions W = C sin 2q7Z 
where C is a constant representing the initial magnitude of the disturbance 


and p, q are integers (¢ ~ 0): thus their general solution can be written 


w—c Mz. (19 


sin 
where 1M = mz and m is an appropriate integer. Although there is little 
doubt that the pattern of convection will be different when N is not zero 
or when the motion grows in a layer which is not uniformly stratified, it 
seems reasonable to assume that these differences will be relatively small. 
There will still be a similar tesselation of the flow in horizontal planes, and 
cells will continue to occupy the whole region between the bounding surfaces, 
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for it seems that any other arrangement would increase the constraint on 
the motion. Therefore, this solution (19) may be used as the basis for an 


pproximate investigation of the dependence of N on the dimensionless 


parameters A, 6, and o; if it is substituted in equation (18), then 


(K2+ M2)N2+-(o!+o-!)(K2+ M2)2N-+ (K2+ M284 K2AB = 0, (20) 


cos” 


where B 2 b(Z) . ,UZdZ 


sin 


s the ‘entire profile number’. 


A simple check can be applied to the present treatment at this point. 


For marginal stability NV = 0, and hence from equation (20) 
> l r9 9\<6 
AB — ...(K?4+ M?)3. (21) 
kK? 
In particular for a linear temperature profile B —1, and in this case 


21) reduces to the criterion for neutral stability given by Pellew and 
Southwell. However, for other values of B representing different profiles 
the critical Rayleigh number will be decreased or increased according as 
B Ol l. 


5. The effect of non-uniform stratification 
It at time ¢ = 0 the lower surface of the layer is suddenly raised in 
temperature from 7; to 7), the temperature profile at any subsequent time 


can be expressed in the form 


o—Z; : 2 — sinrxz(1+Z oer zict 
at ($—Z) S (2 “exp -= ——]. 
Al ° TT —d r h2 
r=1 
The corresponding value of the entire profile number is B = —1. Thus, 


to the present order of approximation, the stability of the layer is not 


iffected by this method of initial heating, although inaccuracy in this result 
Svaeecis ' £008 seu 
will arise from any considerable departure of W from the form C . MZ. 
sin 


By considering a rather wider range of physically possible profile func- 


tions, it may be shown that the entire profile number 


) 4 IT l m—1 (r+)m 1T 
. 2 ( cos- z P ( u 
B 7 = . wAWZdZ li 4 cos u du 
AT J} dZ sin? TAT aw du 
r m rm 
lies within the range —1 B < 0 for all functions 7(Z) which are such 


that (i) dT'/dZ 0 and (ii) d?7'/dZ? > 0 for L<Z: 1. These con- 
ditions are equivalent to stipulating that within the layer conduction is 
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directed from the hot to the cold surface at all intermediate points, and 
that there are no internal heat sources. Thus it appears that the critical 
Rayleigh number for marginal stability will certainly not be appreciably 
reduced in free layers which are non-uniformly stratified, and it may actually 
be increased. This result is of some interest in problems of thermal stability 
because it has been suggested (Sutton (8)) that such non-linear temperature 
profiles may be specially effective in the initiation of instabilities. 

One of the predictions of the Rayleigh theory for marginal stability in 
fluids which are uniformly and unstably stratified is that very large tem- 
perature gradients will be needed to set up cellular convection in thin layers. 
While previous experimental work had shown excellent agreement with the 
general results of the Rayleigh analysis, there had been no special attention 
to thin layers until some experiments by Chandra (9). He verified the 
original results for Bénard cells, but he also found evidence for a new kind 
of convective instability which could grow in thin layers between rigid 
surfaces at Rayleigh numbers which, under some conditions, were far below 
the usual critical values. It was this new type of instability (described by 
Chandra as ‘columnar’ ) which Sutton suggested might grow from a thermal 
boundary layer near the lower surface in the fluid. 

Although the present analysis has been developed for the freely bounded 
layer, the two types of system are sufficiently close physically for the results 
in each to be of similar nature. If the effects of both viscosity and thermal 
conductivity are taken into account, it has been shown that reductions in 
the critical Rayleigh number for the free layer may be expected to be rela- 
tively small: it may reasonably be inferred that the equivalent reductions 
for the layer with rigid boundaries will not be as large as those reported by 
Chandra (ranging up to a reduction of > 90 per cent. for a layer of air 
c. 5mm thick).+ Thus the effect of a thermal boundary layer on the equili- 
brium of stratified layers appears unlikely to represent the true source of 
the reported ‘columnar’ instability; if this is a genuine instability of the 
layer some new mechanism will probably have to be found. 


6. The unstably stratified layer 
The results of this section may conveniently be discussed in terms of a 
modified Rayleigh number A* = AB, and for the unstably stratified layer 
A* <0. The behaviour of N in terms of the non-dimensional parameters 
A* and o can now be found from the solution to equation (20), 
N = (K2+¥?)! (oi +o-4)4 (oto Bo :|*}, 
(K2+M?)3| | 
+ Moreover, there appears to be no reason for anticipating Chandra’s observation that the 
depression of the critical Rayleigh number increased with decreasing thickness of the layer.‘ 
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where the positive sign has been chosen since this corresponds with the 
instability. For given values of K, A*, and o the largest value of N corre- 


sponds with the first harmonic mode, for which m 1 (i.e. M = zm), and 


4 T T T 








69 














4 1 1 
Fic. 1. The initial rates of growth (N,) for disturbances of wave- 
number (K,) in unstably stratified layers, plotted in dimensionless 
form for a range of values of the Rayleigh number (A,) and for 


o l 


since this analysis applies only to the onset of convection other modes need 


not be considered. If relation (22) is reduced by the substitution 


M TT, K, K |x, A, —A*/74, N, = Nis", 
it takes the form 
| a viene SKIN 
N, s(1-+ AY), (ao! +o-!)+ | (ot—o*)? 4 3 ma \ (23) 


The variation of the rate of growth coefficient (N, = 4?n/7?(«v)!) with 
the wave-number (K, = hk/z) is shown in non-dimensional form in Fig. 1 
for a variety of values of the modified Rayleigh number (A, = —A*/7‘*) 
and for Prandtl number o 1. Although the present linearized theory 
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applies only to infinitesimal disturbances, the plotted curves do serve to 
illustrate the earliest stages of the growth or decay of disturbance modes, 
There is a critical Rayleigh number A,, = 2? below which all wave-numbers 
will be damped out, while for A, > A,, there is in each case a band of waye- 
numbers for which disturbances will be amplified. In such cases of instabi- 
lity one mode will grow more rapidly than the others, and it can be seen that 








-4 








“2 


Fic, 2. The dependence of the rate of growth factor (N,) on the wave number (K, 





for various values of the Prandtl number (c) in the case of marginal stability (is 
A, = #7) for an unstably stratified layer 


the wave-number of this dominant mode increases slight ly with increasing 
Rayleigh number. Thus disturbances with larger wave-numbers will begin 
to grow in those layers of fluid in which the effects of diffusion are small. 
The broken line which joins the maxima of the set of curves in Fig. 1 has 


the equation . oe 
ms 1+ Aj+ 2K}. 


These results are in broad agreement with observations, although variations 
in the dominant wave-number have not been (and are unlikely to be 
observed because of experimental difficulties. They may be contrasted 
qualitatively with the equivalent results of Hide, according to which all 
disturbance modes are unstable. 

The symmetrical way in which viscosity and thermal conductivity affect 
the results is further illustrated in Fig. 2 by a set of curves showing the 
variation of the rate of growth factor (N,) with the wave number (Kj) 
for several values of c in the case of marginal stability (A, = 27). It can be 
seen that the critical wave number (K,) for marginal stability is unaffected 
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w changes in the Prandtl number. However, this is not true generally, 
nd for other values of A, the dominant wave number does depend on o 


osome extent: this effect is not due to the predominance of either viscosity 


r thermal conductivity but only to their lack of balance, since the effect 


sunchanged by substituting 1/o for o. 


7, The stably stratified layer 

Attention has been concentrated so far on the equilibrium of unstably 
stratified layers where the neglect of thermal conduction in the fluid is likely 
to have the most serious effects. In order to complete the qualitative 
mparison with the treatment given by Hide (6) a brief discussion will 
w be given for the stability of stably stratified layers. 

In this case A* 0, and if A, is now written for + A*/z* the corre- 


sponding solution for N, is 


i 10 4KEA, 13 
N, L(] K2)! (o o~*)+ | (a'—o-*)° 1 | (24) 
: | (1+ A?)3] J 
ence, N, will be real or complex according as 
F(K (oO oO )? 4K} 0 
A, (1+ 4G)8 
[he equation F(K,) = 0 has either two real positive roots (,A,, .K,) or 


none at all. In the former case if A, lies outside the range (,K,, .,) it 
follows that F(A,) > 0 and that JN, is real and negative: there will conse- 
juently be aperiodic damping of the corresponding disturbance modes. If 
K, < K, < .K, then N, is complex with negative real part: in this case 
there will be periodic decay of these disturbance modes. Finally, if 
F(K,) = 0 has no real positive roots, N, will be real and negative for all 
values of A,, and hence all disturbance modes will be damped aperiodically ; 
this will be the case when A, < 4(a!—o~!)*?. These results differ from 
those of Hide in numerical detail only; such a difference is not significant 
is the present results are for a slightly different system and are to be 


regarded as of qualitative rather than quantitative value. 


APPENDIX 


A COMPARISON OF THE RATE OF GROWTH OF THE 
DISTURBANCES AND THE RATE OF DECAY OF THE 
TEMPERATURE PROFILE 
If the form of the disturbance is represented by w C's, (mare h) exp(nt) it 
follows that the time in which its amplitude (|w]|) grows by an order of magnitude 


Ss approximately t, 2-3)n 2-3h? (772(«v)*N, ) seconds, where N, (> 0 here) is the 


limensionless growth factor which has been plotted in the figures. 
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In order to find the change that will take place during this time in the non-linegy 
profiles it is necessary to consider special profile forms. Two of these will be sufficient 
for the present purposes: 

(i) The profile considered already, and caused in a layer which is initially at 
uniform temperature 7 when the lower surface is suddenly heated to Ty (> 7) and 
maintained at that temperature. If the time is measured from this instant, thy 
profile at any subsequent time is 


x 


T —T, 2 N¢ sinra(44+Z Trt 

eat ee a Sinn! Bad mialidlen f 

iA r s 1 hud r h? 
r=1 


(ii) The more realistic profile, which is caused in a layer initially at uniform ten 
perature T,, when, at ¢ = 0, steady heating of the lower surface is suddenly started 
at the constant rate k(7,— 7',)/h units of heat per unit area per second. The ultimat 
temperature of the lower surface will be 7), and the profile is 


T — T, , 8 (—1) P 
(}—Z) S sin $7(2r+-1)(3 Z)exp| 
— 


qr?(2r+-1 yin 
™%—T, ° 7 (2r+ 1)? = 4h? 
7 


In each of these cases the ultimate profile of temperature across the layer is 
T’—T, = (T—T,\4—2); 
for case (i) a simple and approximate estimate of the non-linearity of the profile is 
given by the value of (7’— T)/(T,— 7) at Z 0. In case (ii), however, the tempera- 
ture of the lower surface, ;7) say, increases gradually to its final value 7); thus th 
corresponding measure of non-linearity must be taken as (7”— T')/(;7,— T,) at Z = 0 
The times (¢,,) in which ( 7” Z) (Z. T,) and (T”’ T) G7. T,,) are reduced to one-half 
of their initial values are, rather approximately, h?/7?« and 4h2/7*x, respectively. 
These times, t,, serve as a rough time scale for the decay of the non-linear tempera- 
ture profile, while ¢, serves as a time scale for the initial growth of the disturbances, 
The ratio ¢,/t, then provides an estimate of the extent to which disturbances can 
grow before the temperature profile has changed significantly towards the linea 
profile corresponding with the steady state. Substituting the results obtained, 
t,/tg = 0404N, or 1-604N,, respectively, so that a value can be given to this rati 
for each dimensionless wave number (K,) and Rayleigh number (A,) from the results 
which are represented in Fig. 1. It can be seen that only the more unstable modes 
will grow markedly in the time available, and that for A,—?, small, none of th 
unstable modes will have time to be affected greatly by the special features of the 
temperature distribution. Moreover, these estimates will allow for rather more growth 
than is likely to occur in practice, as the exponential amplification continues only 
while the disturbance is still small and subsequently a steady pattern of convection 
is realized—in this case Bénard cells. It may be noted that Thomas (10) has alread) 
pointed out in a discussion of Bénard cells that strongly non-linear temperatur 
profiles cannot persist in layers of real fluid under normal experimental conditions. 
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THE MOVING AEROFOIL IN SHEAR FLOW iN THE 
NEIGHBOURHOOD OF A PLANE BOUNDARY 


By D. E. EDMUNDS (130 Stanwell Road, Ashford, Middlesex) 
[Received 25 October 1956] 


SUMMARY 
In a previous paper (1) we have discussed the problem of finding the forces acting 
on a two-dimensional aerofoil moving in inviscid, incompressible fluid bounded by 
a plane wall, the fluid being at rest at infinity. Here the theory is extended to thy 
case where the undisturbed fluid motion is that of linear shear flow parallel to th 
wall. Expressions are obtained in series form for the forces and couple acting on; 


fairly general type of aerofoil, and the flate plate aerofoil is considered in more detail, 


1. Introduction 

Coomss (2) has dealt with the question of a stationary aerofoil in linear 
shear flow in the neighbourhood of a plane boundary. He uses the method 
of Green (3) which was developed for potential flow, and expresses the 
results in series form. James (4, 5) has discussed both the problem of the 
stationary aerofoil and that of the moving aerofoil in shear flow, but in 
an unbounded region. 

In this paper we consider an aerofoil moving near a plane boundary in 
fluid which is in linear shear flow parallel to the boundary. It is supposed 
that the fluid is inviscid and incompressible, and that the aerofoil is two- 
dimensional. This problem does not seem to have been discussed previously, 
and it is found that the extension of Green’s method used in (1) may be 
employed. The expressions obtained for the forces on a general aerofoil are 
complicated, and so detailed consideration is given only to the simple case 
of a flat plate aerofoil. 


2. Construction of the complex function Q 

Consider the two-dimensional motion due to a cylinder moving near an 
infinite plane boundary in fluid whose undisturbed motion is that of linear 
shear flow parallel to the boundary. This motion is considered in a plane 
defined by a normal cross-section of the cylinder. 

Define two sets of moving rectangular cartesian axes, the (x, y)-axes 
fixed in the aerofoil and moving with it in all respects, and the (2’, y’)-axes 
whose origin is fixed inside the aerofoil section but which do not rotate 
with the aerofoil and so remain at all times parallel to a set of axes fixed 
in space. The axes fixed in space, the (a, yy)-axes, are so chosen that thei 
real axis and the boundary wall are parallel. Thus the plane boundary has 


[|Quart. Journ. Mech and Applied. Math., Vol. X, Pt. 4 (1957)] 
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MOVING 


equation y’ -b at any instant, and the force components parallel to 
the 2’- and y’-axes may be written as X’ and Y’ respectively. Making 
the (x, y)- and (a#’,y’)-axes have a common origin O inside the boundary 


urve C of the aerofoil, the relation z’ — ze“? holds at all times. Here 6@ is 


the angle between the (x, y)- and (x’, y’)- axes at any time, and z = x+y, 


a’ -+-iy’ in the usual complex variable notation. The various coordinate 


svstems are shown in Fig. 1. 


Yo 


























—_— 
b 
y' a) 
H 
| in 
Fic. 1. The coordinate systems 


Referred to the (xp, ¥)-axes, the undisturbed shear flow considered is 
u(l-+ykl), V=o0, (2.1) 
where U’, V are respectively the a, and y, velocity components of the flow, 
,and k are constants, and / is a typical length in the aerofoil (e.g. the 
hord). 

With respect to a set of fixed axes which coincides with the instantaneous 


position of the (a’, y’)-axes, this flow has 2’ and y’ components 


uo(F+-ky'l-1), 0 2.2) 


respectively. Here F |+-Bkl-1, (a, 8) being the coordinates of O referred 


to the (2, y,)-axes. 


\s usual, a stream function y% may be defined in terms of the absolute 
elocity components of the fluid parallel to the 2’- and y’-axes by 
ods , o ‘ 
- Fee. (2.3) 


oy 6x’ 


™ 


the Helmholtz equation for two-dimensional motion is 
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¢ being the vorticity vector, and 7 the time. In terms of % this is 
C o op C¢ " 
. te ‘ |) Vb = Q, (2,4) 
da’ Cy’ = ey’ ex 
Choosing the undisturbed stream function %, before the introduction of 
the cylinder to give 4% = 0 on the plane boundary, we have 
hy u(y’ +b) F+tk(y—byl-}} 


oe le k 
so that Vib : bo, obo = — We, (2.5 


da’'2 oy’? =f 
The obvious solution to (2.5) and (2.4) is that the vorticity remains constant 
and equal to u,k//, i.e. % must satisfy 
Vb = uy k/l. (2.6) 


The uniqueness of this solution follows from the assumption that the 
presence of the cylinder produces no change in the total vorticity of the 
fluid. 

Following James (4) we introduce a stream function ¢’ defined by 


b= b+’, (2.7 


where Vb’ = 0. (2.8) 


Hence y’ defines an irrotational disturbance motion, and a complex poten- 
tial funetion Q’(z’) may be defined which has ~’ as its imaginary part. 

Define also a complex function 2 of which the complete stream function 
is the imaginary part. The real part of Q is not a velocity potential since 
the flow is not a potential flow. 

If the aerofoil has a translatory motion defined by velocity components 
u’, v’ along the 2’-, y’-axes respectively, and has an angular velocity 
about O, then ’ satisfies the conditions: 

(i) Q’ is a function of z’ but not of 2’ (a bar placed above a quantity 
will be used to denote the complex conjugate of that quantity throughout 
this paper). 

(ii) On the aerofoil boundary 


. , 7, y a 9 '-r —s 7 t-/ 
im) +uU, F2’— ity U(2’'2?—2'2')+-W'2’ —hiwz2'2 | constant, 
where w’ wu’ +iv’ 
(iii) im(Q’) = constant on y’ aah, 


(iv) Q’ reduces to zero at infinity, where the effect of the aerofoil is 


imperceptible 
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The function 


, =. (A A 
()’ iu, A,{logz’ — log(z’ + 27b)}+-uy S a ee (2.9) 
— ‘ 


satisfies conditions (i), (iii), (iv), and includes a circulation K = 27A up 


ound the aerofoil. 
We define the aerofoil considered here by the transformation 


z=t" > af, (2.10) 
n—0 
vhich maps the boundary of the unit circle |t | in the t-plane on to the 


joundary of the aerofoil in the z-plane. The region outside the aerofoil C 
transforms conformally to the region inside the unit circle y in the t-plane. 

As the method of solution is now similar to that of Green (3), only the 
essentials will be given. 


In terms of ¢, (2.9) becomes 


iuy Aglogt+uy > (C,,t"+D, t-"), (2.11) 
n=1 


where the C,, and D,, are defined in the Appendix (section 1). Hence the 
oundary condition (ii) leads to the relations 


uC (i Up Fy wa, 1 


Uy D,,+iwb, + ikuy(?y,—2b,)l-  (n > 2) 





Ug Us (w+ Uo F )e "Ay (2 12) 
Uy Do+twb,+- fikug(*y.4 29 _»° 2b,)l . 
Up ( " (iw Up F )e “A, 
Uy D, iwb, T (Ww Ug F ye %G, T Likug(*y, 1 2 1 — 26,)l- 
where "y, and 6, are given in the Appendix (section 1). 
The wall y’ b will have equation y, = H referred to the (xp, yp)-axes 
fixed in space, H being a constant satisfying 8 = b+ H, so that 
F = 1+$kl-1 = G+6kl-, G = 14+ Akl-. (2.13) 
Thus the boundary conditions (2.12) suggest expressing A, in the form 
A, > "A,b- (n> 0), (2.14) 
8 1 


ind the C,, and D,, can then be written as power series in b (Appendix, 
section 1), Substituting these expansions in (2.12) and equating coefficients 
of corresponding powers of b give’ a set of equations defining the coefficients 
A 


this determines Q’, and hence Q, except for the arbitrary circulation. 


n > 1)in terms of the circulation coefficients °A, (Appendix, section 2). 
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. Th 
3. The forces and couple on the general aerofoil 
The force components Y’, X’ and the couple I about the origin O are 
given in terms of the pressure p by the integrals 
Y’+24X’ = p dz’, I’ = re | pz dz’, (3.1) with 
. o 
taken round the boundary curve C of the aerofoil. The ordinary pressure 
equation cannot be used since Q is not a potential function, but p is deter- 
mined by using James’s method (5), which is given in some detail. so th 
He returns to the equation of motion of the fluid, i.e. to p 
D l 
a -—Vp, (3.2 p 
Dr 
p being the fluid density, and q the vector velocity of the fluid, referred to 
axes fixed in space and coincident with the instantaneous position of the 
(x’, y’)-axes. Writing (3.2) in component form and using (2.3), gives 
ou oA oV OA 
or Oa”” or ey’ | dots 
‘ (3.3 U 
9 9 79 ) 
where A= ¢V*2s—4(U2+ V2) | 
p 
Equation (3.3) is equivalent to 
fj ou ’ ou whic 
C Ou a . OW a 
Eri 0, - (arid) — 0 (3.4) | and 
Cz OT CZ OT . ge 
indic 
which has obvious solutions 
y'4 
. Os ?  : — oe : 
=i pe), ra=—i%s pee, (3.5) | pi 


C7 CT 
f(z’), f.(2’) being two arbitrary functions of 2’ and Z’ respectively. Hence 
. Os — Se 
21 = fo(2)—fi(2’). (3.6 
: 


C 


Now ¢% = +’, where éy,/er = 0, and yy is the only term in ys which con- 
tains products of z’ and 2’, and which therefore cannot be split into two 
functions of z’ and 2’ as in (3.6). Thus a solution of (3.6) is possible, and 
9 Op o 0’ _O’ ~ ie 
21 = 5 (Q’—) = f,(2')—fi(2’). 
OT OT 


Since 2’ is a potential function, 


ee 0Q’ 
fiz’) = fans , 


which gives 
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The pressure equation is then 
) 7 i, wo , Aeon |, Oe 
I bb Vd 1(( at | J “)+ ( 4 . 
p : ale ; 


vith V2e ujk/l. It is easily seen that 





p h ’ 1/6Q  é6Q\/eQ = eQ\ /dz dz\- 
vot aa — sala a 
i) / 2\ ot ct of ct dt dt 
ae ee 0’ [dz\- = . ., OG)’ (dz\-1 
L(w’ et iwz | 4 (w'e" —iwz) — | -| (3.8) 
ot dt “ ct dt 
l C ( c as = CO \ P , 
WwW i " u , u an. WwW (Q2 ()’) 
2\ o6 ob Ow OW Cw 
ts denoting differentiation with respect to 7. 
Using (2.11), (3.1), and the relation 
Q = QO! 4-u, Fz'—Jikug(2'2—2'2’-2 (3.9) 


vhich follows by definition of the complete complex function Q, the forces 
nd couple may be calculated. The various integrals are evaluated as 
ndicated in (1), and it is found that 
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o 00 > ) 
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Up W i A,a, D, a, > na, 1Cnyy—- 
n=1 
iw, tA (bo lo) > nC, l- -D,, l, T C, b, D,, b,,)| ass 
? l 
¢ i. c Cc 7 C 
Up| WwW u Uu o-— v 
ou cu Cu CW ob 
121A, > b ¥ (C_,+D,) ¥ (r—1)a,,.,4,)4 
n 1 n=1 r=0 
ae . - 
> (¢ n't D,,)\ Zz (m wg 1)a,, ra, i 
n=1 r=0 
where 
c C,,+ Fea, +1; D; D,+ Feay, D;, D, (n> 1) 


nd the coefficients c,,, 7,,, d,, T,, Rns Ras bas Us fins fins Mn» me, are defined in 


the Appendix, section 3. 
The form of these expressions suggests writing 


1 ¢X.)b-8, 4 % £0, (3.12) 


2 s=—2 


Y’+iX > @, 
ind these expansions can be obtained for any particular aerofoil. 

To complete the solution a definite value for the circulation must be 
obtained. For an aerofoil with a sharp trailing edge 7’, this can be done 
by assuming the Joukowski condition to hold, so that the fluid leaves 7' 
smoothly with finite velocity. The velocity at 7' due to the term yy in & 
is clearly finite, and the transformation (2.10) can always be chosen so that 


T corresponds to ¢ |. Hence the condition requires that @éQ’/et = 0 
it t |. Proceeding as in (1), the circulation coefficients are found to be 
ri - (3) } ] 
A, ik(d, ¢ a,e”)i- 
1, Uy slUy” 1 (Lp € Ay ad 
if(w! +p Bea — (i +g @)eia}-+-w > (—1)"nb, 
n=] 
| j = 4= ™ os 4 
Ugh \A (Ao ay je“ "+ Ay(Qo A )e nae > ( 1)"nb,,|1 ,, (3.13) 
1 
1 1 ” 
( { Us } > > ( 1) ly 0 - ny : m Xe 4 n'y) 
1” ? 
1 7 1 wii 
NN + N ] ms n s-—r ma, s-—r 
é cad a 1 — \ l) nu / i "Og rT A, m ‘ a Os r a - 





4. Discussion of the solution 
The work of this paper differs from that of (1) in the appearance of 


terms involving 6*! and b*? in the expressions (3.12) for the forces and 
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couple. This is because at infinite distances from the aerofoil % reduces 
to w%, which contains terms of this sort. On putting uw, = 0 these terms 
in 6+! and 6+? vanish, and the problem reduces to that treated in (1). 
By replacing H by (8—b) wherever it occurs, it may be verified that 
(3.12) can be expressed in a form which contains terms in f and 8 in 


9 


place of those in }+*, b+1, H, and H?, i.e. as 


Y’4iX’ =D (V%4iX7o, T= >. (4.1) 
s=0 s=0 


The assumption is of course made that these rearrangements of terms are 
valid. 

The expansions (4.1) are often more convenient than the alternative 
expressions (3.12), and will be used here. In particular, they give a clearer 
picture of what really happens in the limiting case of the aerofoil in un- 
bounded shear flow. Using (3.12) the limit is deduced by allowing 6 and H 
to become infinite in such a manner that (b-- H) remains finite and equal to . 
(4.1) provides the limit merely by letting 6 become infinite, when the 
forces and couple are Y¢, Xo, and I5 respectively. 


5. The flat plate aerofoil 
The plate is defined by the transformations 
z = 4U(t-!+4), z’ = ze, (5.1) 
where / is the length of the plate, and @ the angle between the plate and 
the wallat any time. Expanding the forces and couple in the form (4.1) gives, 
after some reduction: 
Yo/7p = Uu'+-uy F){(u'+u, F)sin 6—v’ cos 6}— 
— thluy(2v' —lw cos 0)sin 20— 11?(%' sin 6—#’ cos 0+- $1) cos 6 
Y5/7p —11*{(u'+-u,F)sin @—v’ cos 6} 
x {2(u’+u,F)sin 0—v’ cos 0+ wl(1-+sin? 6)}+ 
+- $y kl*{(w’ +-uy F)sin 6— }wl cos*@}sin?6 4 
-+-(41)8{2(a’ sin 0@—#’ cos 0)+-4al'sin 20 
Y3/7p = 7 2{(u'+u, F)sin 0—v’ cos 6} > (5.2) 
x {(u’ +u, F)(1+3 sin?@)— 2v’ sin 26}- 
gwl"{(u’ +-U_ F)(1—9 sin?6+-4 sin*6) + 
+ $v’(9—4 sin?6)sin 20} + 4w?(4/)5(4—3 sin?@)sin 04 
+ 4(ugk?/1?)(41)> sin®@— (uy k/1)(41)! sin 0f16(u’ +-u, F)sin?é+- 
+ 8v’(1—3 sin?@)cos @—wl(2—9 sin?0+-5 sin‘6)' — 


— 2(41)*{(w’ sin @—#’ cos 8)(3— 2 sin?6) + 3a cos?6}cos 6 
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‘1 yu, F)sin @—v’ cos 6}+- 


Lu, kl(2v’ —wl cos @)sin?6 + 217{2(a’ sin 6—b 


11?{(w’ +-u, F)sin 6—v’ cos 0—} ugk sin6} 
(2) pol cos @)sin 6 

3, 3(a’ sin @—0’ 
F’)sin 0 


S(w’ tu v’ cos 6\cos?6 


uy F)(1 


6w?(4/)> sin 6 sin 20+-4w(41)4{(u' 4 


Sr 
+-v (oO 


Luh 1)(41)4{ 80 cos @—wl(2 
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w( 11)? cos @ sin 20 


2(41)4(w’ sin 6 v cos 6)(3 2 sin?@)sin 6 


126 , ry 


u, £')cos 6 v’ sin 6} 


{(u’ +-u, F)sin 6 
v'(3 


ds Uy kl? sin Of (w’ Uo F’)sin 26 


wl eos 6 


ig/°(w’ sin 6—0’ cos 6)—6a(41)* 


Tp i, 1?{(u'+-u, F')sin @—v’ cos 01 


{(w’ 


sity KU sin?6(2v’ —fal cos 6+-1u,k sin 
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34 sin76)cos 6+ 


0 
w(41)*{(u’ +u,y F)(5 

+y'(37— 
uy F)(1 
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(u_ k/l)(41)° sin Of3(u' 


v’(1) 32 sin*@ 


6w*( 1])6 sin 26 


(41)°{(a’ sin 6—@’ cos 6)(5—2 sin20)-4 


th # 1+Bkl-1. The restrictions l/b < 1, k - 


‘pansions will be more likely to be convergent. 


SHEAR FLOW 


‘ cos @)+-la}sin 6 


4 sin20)sin #}sin 6+- 


5 sin?6)'sin @ sin 20-+- 


1 °c sn YAR 
lu, ksin @sin 26} 


+- Ue F’)sin 26 


-4sin?0)sin 26 
-~5 sin?6)cos @— 


hu, k sin @ sin 46} +- 





cos 6-+- 11a)sin*é 


4sin?@)cos6+- 





v’ cos 8@— wl} 


2 sin?6) + 


v’ cos 26}4 
) S 
6)+-1al}sin 6 


@sin 2 


10 sin?@)sin 6}— 


34 sin?@)sin 6}+- 





1 @1(3—2 sin?6)} 


1 are made so that the 
Further terms of these 


‘\pansions can be found if desired, but the reduction tends to become 


volved and tedious. 
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On putting wu, = 0 these results reduce to those obtained in (1) for the 
plate moving in fluid at rest at infinity. Agreement with the force expres. 


sions given by Coombs is reached on writing wu’ = 0, v’ 0, w 0, F =] 
> e S 


in (5.2) and (5.3). He does not give the couple, but (5.4) shows that this 


is given by 

— >, = (14+ $hsind— 4k? sin?6) — 3(1/b)(1— 74k? sin?@)sin 0-4 

7 pl*us sin 26 ? 5 

+-g5(2/b)*{(1+- 10 sin?@) + 3k(1—4 sin?@)sin 6 . 
—},k? sin? cos 20} + O(b-%) 

For small 6, 


D/P) = 1—4}(/b)(1—fsksin 0+ 1k sin?6)sin 64 
+3 (1/b)*{(1+- 10 sin?) + 1k(1— 32 sin?6)sin @— 1k?(1—38 sin?6)sin?6\+.,, 


neglecting powers of (k sin @) higher than the second. Hence I'/I} increases 


with k, for positive 6., Thus for @ 10°, 
r i‘: 1—0-1158(1/b)+0-0519(1/b)?-+-..., 
T/T = 1—0-0868(1/b)+-0-0407(1/b)2- 
P/T) = 1—0-0695(1/b) +-0-0340(1/b)?- 
for k —0-5, 0, and 0-5, respectively. 

In this case the shear flow is associated with the aerofoil in that th 
condition 8 = 0 is imposed. Removing this restriction the results tak 
the more general form 

Y’/Y9 = 1—4(U/b){1—(k/4 F)sin 6'sin 6+ 
+- +h. (1/b)?{(1+-3 sin?6) — (k/ F)sin36+- (k?/16.F?)sin46\ 4 
Yo = mplu? F?sin 0 
A = 
aan 5p = (F2+ 4 Fkesin 6— yh? sin%®) 
$(1/b)( F?— +k? sin?6)sin 6-+-3,(1/b)?{ F?2(1+- 10 sin2@)4 


+? Fk(1—4sin?@)sin 6—+,k? sin?6 cos 20}... 





By giving f different values, it is now possible to determine the variatio! 
in Y’ and [ as the plate is placed at different points in the field of flow. 


teturning to the general expressions (5.2)—(5.4) and writing v’ = 0, =! 


u’ = —u,F, gives 
= stam pus k®1(1/b)? sin®@+-O(b-3) | 
x =D ‘ie 


D/ap = —jhguj k°l* sin26 sin 26{1 —4(1/b)sin 9+-3,(1/b)? cos 20-+...} } 
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Here the plate is in motion parallel to the wall with velocity exactly equa 
0 the undisturbed ctream velocity at the centre of the plate. In an 


nbounded fluid the forces acting would reduce to zero, but the presence 
fthe wall means that the small force Y’ acts on the plate in a direction 
erpendicular to the wall. The shear nature of the undisturbed flow 
usures that a couple acts even in the limiting case when 6 is infinite, and 
» notice that for small positive 6, the magnitude of the couple first 


lecreases and then tends to increase as //b increases. 
Finally, consider the case when the plate executes translational oscilla- 
tions in the direction of the y-axis. This is obtained on setting 


’ ‘ . 
uM Q. W 0. v VU, SIN NT, 


and » being constant. It follows that 


v v : 
S = Q 0 « > 
b h(1 “7-cos nz, p p(1- = cos nz} 
. Po j 
shere the constants / and £, are the mean values of b and £ respectively, 


nd 8, = A+H. Writing F, 1+k8,1-' the forces and couple become 





22 k2y2 \ (1 
—= ui sin 6| | Fe i F?—1kF sin 64 : \ sin 6+- 
oat : | as 5\ iil 21?n2/\h 
1 // = ? ’ kv2 ‘ ‘ 
F2(1+-3 sin? °- (k+-16F, sin 6+-k sin?) + 
amar ol me) Open? | datas, 


js (k sin 6—16F))k sin | = 


- ya(-)(1 - sin 6 }cos*6 


3/1 
— U,V, cos @sin nt| (f+ $k sin @) — | | sin 6+ 
. | 4\h} ° 


1 /l\?(,, 3kv2 . ,) 
(7) | F\(1+-7 sin?@)+- $k(1—3 sin?@)sin 6 + mp in| sie 


red F, sin 6) — 


> 


7 
+5 nvgeosnr| (cost 


ris) a , — ve en. 
cos? — (8 F,—k sin @)sin @}sin 6 ——° cos?6+- 
lh} | Pn? hehe 


(1\2 —_ ss 
(3—2sin*6@)cos*6é 


| 
30\h 


2 sin 6 =e : — spends 3k2v2 . 
it 14 * sin 0+kF(1+-2 sin?é)— 2k? sin*@ +-—; sin a\| 4. 
22 


2n*h? | 
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52 
_ Uo % 


; cos 6 sin 2nr| 4k 3(; isin 6+ (;) (Fosin? + ¢k(1+-7 sin*6)} | + 
v oy 


on 


9 


+ 4v% cos 2nt 





| Pn? 


3 2 
Yo {2In(cos%9 an up 
~ 32n2h?| [?n? 


a , 1\2) 
— 1y*cos 0/1 - 

a. (a) er rata) |— 
i 


: . . ; | 
~ wg rosin Osin nr] (Fy isin 6) A) 


t 


l me 1, 
— — Nv, sin 20 cos nz! I ( }sin 6+ | 
8 A\h 3: 


j 1/1 
k. | b sin O-- 
; i) in é 


9 
UnVy . ° 
ak ° ° sin @sin 2nr 





L6\A 


: ijt\. , 1\? 
— Jnf.cos 8 cos 27} i(,)sno+ ali) , 


7% 
+ sin 4(/n cos 6 cos 3n7— kuy sin 6 sin 3n7) 4 
i6n2 h* 
and 
Poy lw: kv 
—. = }u?sin 20] | F2+-3h(F,—}ksin 6)sin 6 + 
m7 pl? | 


L ft\*{,. , oe . fe 
+ sa(i) |(1-+10sin%6)( 3+ TE a) 


« 7 . 9 ° 16k al 2 
-3kF,\(1—4sin?6)sin @ + an F, ve sin 6— 





OM ind 4 (cos? 26 - neha 
ay* ?n? h 


ie 16 
— “(;) sin 6/5 cos0 0 LF s in 0 —~ 2 


sin‘ *0)cos 3n7— 3u,k sin 26 sin a + O(h-) 


;) (3 2 sin29) +. — 


] .. ; ‘ 
+ (;) {4 sin @+-k(1—3 sin*6)}| - 


7 = 
a 


I l 12 9°. 9 . k2y2 , 
— ve oO is k?sin?0 + Py) §+ 
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}cos 20... az) 10 sin?@)sin a _ 





XN 


l 
' 12 cos 6) sin @ | | 
| 4 


1 ; 
+ Uo U9 


sin nz| { F, cos 20+-14(3—2 sin?@)sin 6}— 


| (Fi, sin 36+-3k sin®)(_ 4 salj) (0 +18 sin?@— 20 sin*@) + 


nail 
. 


—!— sin 36 | +- 
l6n7h? 


1k(11—32 sin?6+ 24 sin4)sin 6} + 


l 
=e NV, COS 6 COS NT 
16 


F 2hu? 


| l?n? 


(4F,-+-k sin 6)sin D - 


l l6us l 1 /l\2 gy2 
— i °kF sin 6 sin 6+ —{ 5—2sin?6) —_ —°_cos 26— 
| oO i) sala) irate 


1 4n*h? 
“oO sin 6 2( F3 on O\ sin 0+ }kF(1+ 10 sin20) + 
n*h* i \ 4]? 


+g k?(3— 16 sin?6)sin a + 


Sin 


M0" 0 sin 2nz| kos 26 bi( sin 36 ! 
h 


l > ) (Fosin 360+2k(1 + 22sin*#—20 int 4 
t 


| hud)... Al ku? . 
1a? cos 8 cos 2nt | ie “*)sin 6+ | cos 26 — - “0 sin?) — 
' ANA 


> ) 9 9 6 9 
_ [2n2 12n? 


L (l\? ku ‘ 
0 sin°@)sin 6 + 9 {39F sin 6+k(1+10sin26)'s ' 
sz) (11—10sim*@)sin ¢ 39n2h2? F, sin 6+-k(1+-10 sin26)}sin a af. 


a] k2 2 
t v"Up . : . 
+—— 5\In{ cos 8 cos 26 ——* sin @ sin 20 }cos 3n7— 
64n2h?| [22 


—ku, sin 36 sin 3nr\ + O(h-). 


[t is seen that, as in the corresponding case in (1), terms occur in these 
expressions which have frequencies of twice and three times the frequency 
of the motion of the plate; also from X’ a thrust can be obtained on the 
plate. The mean values of the forces and couple over an interval of time 
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27/n follow at once, and denoting these mean values by the left-hand 
subscript m, we have 


ad 


} 9: 719 k2v2 ] -_ ka oe ] ; . 
i usin | (3+ swan) —3(%3 -1kF, sin 6+ Pn }(;,)sin + 





i wala) (730 1-3 sin20) +- =, (k--16F, sin 0--k sin?) -+ : 
+ 45(k sin @—16F,)k sinsol | — _ 
at vali) — sin A) cost O(h-) 
nN = argoos{t— (fino + ECF) + a 
mi = u2sin 24 F24-1k(F,—ksin 8)sin 0+ - 09 | . 
7 pl* | On| 
(5)(¥3 — Pk? sin?6 + ee) 6+ 
i sala) (" { LOsin*#)( 3 rs = ay 4. 8kF,(1—4sin26)sin 04 Pe 
+o Fy v3 sin 06— +k? sin?@ cos 20\| se 


l 1 /i\? — 
v2 cos 6! sin 6 +. s 20 - 7—10sin26)sin 6\ + 
Vp CO in ‘(ie si) (7 sin*@)sin 7 


O(h-*) 
The mean forces and couple can now be determined for any specified 
plate motion of this sort. 
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APPENDIX 


1. The coefficients C, and D, are defined by 


Ch . i°B, Ag T ~ "Bu Ar T "Yn B,, D,, >’ n B,, 
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r r=1 r=n 
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where B, "YX A, b-” 4 > Nee A, b-("+r) 
r 
inti n+r—1\ir7 
and "ty = —>a» a, ( laws (n, 7 1). 








[he ret 


which 


Usu 


" ane 





ft-hand 








.(5.10 


ecified 


advice 


grant 





THE MOVING AEROFOIL IN SHEAR FLOW 


> 


he remaining coefficients °8,,,"B,,"v, are expressible in terms of the a, by the relations 


x x 
> B,t", og 2" logt 3 Bt", (2’)" > "yy t" 
f n=( n -— 
ich follow from (2.10) on omitting terms independent of t. Finally, 
b, Dd an. rn: 
Using the expansion (2.14), and assuming interchangeability of summations, the 
and D becom«e 
| © st+1 
Y 02 * Q \ ~ ~m a 0 8 
( > \y “ . 1 3 ” 1 (5 x > myn - Xo 1 m b 
1 s=0 m=1 
D @e-1 ¢+} 
' " ~ Ma, ™m s— 
— Zz Zz "Yn Xg Ayn b > 
s=1 r=0 m=1 
- s+1 0 e~i 221 
D > 2 > Y x9" b > > > ay ea Aye b 
l a s=n r=n—-1 m=n 
2. The set of equations defining the "A, (n 1) in terms of the °A, is: 
0B,°A_,+18,14_,+e%a, kl“! eG, ki-!, 
8 °A + >'B,"A_,+e%a,,,k0-' = 0 (n> 1), 
of 2) Ly 181A, . 1" Xo. i (w’ Uy G)e "dy 
Ug} " Xo A iwb, (w Uy Ge ‘Ay 
tol "Be" 1, = Bs 1 _ 21 0%9°A J (w Uy G)e as 
1 
iwbs, }iug(*y2 *y . 2b, )kl Ri 
ay? = Fa ae ) ¥9°A_; (w’ + Uy G)e”’ayn., 
1 
iwb, + iug(*y,—2b,)kl-1_ (n 2) 
s+1 s—17 a 
3 > B { > ; ¥,°.4 = = se X< 2 ¢ oe 
1 1 r=0 m=1 
/ e s—1 r+1 
> y Xo" | > NS may Gi 3 r a ‘a (s i. 8 n) 
r=n—-1m=n 
4 1 XA _y (s l,n s+1) 
0 (n s+] 2) 
3. TI fficients ¢,, J», dy, d., hy, hi, ly, Ui, were used in (1), and are given by 


Jn 





n 
> a,¢ 
a 
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Lastly, the following expansions hold on y: 


ee —-1 dz ™ i " s f’ 

zi—)  — t t” 4 t-") 

ai} di n=0 “ n=1 . ; 
heed @ n 

where Ia = FG Egy Sn = Zhi + S46. as 
r=0 r=0 r=1 

anc zz ) —=-l { mt” + m, t "), 
\dt dt n=0 n=1 
DD x n 

where iy: = 2 aR, ., m, = Dak + > a,,h. 
r=0 r=0 r=] 


These define f,, f;,, My, m},.. 
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SUMMARY 


After developing Goodman's solution for radially symmetric modes of stress wave 
pagation in infinite slabs, this paper gives first a detailed analysis of the determina- 
of the various modes. Then, curves of both the non-dimensional phase and group 
given for Poisson’s ratio from zero to half. The form of excitation at 
entre is discussed as well as the limiting cases for thin slabs and surface waves. 
theoretical curves are compared with some recent experimental data from Bell 
phone Laboratories and fairly good agreement was found. Some experimental 
rk at Frankford Arsenal is also discussed from the point of view of the present 
ory 


ry, 0, 2 = cylindrical coordinates. 
U,, Up, U, displacements along r, 6, and z-axes. 

W,, Wp, W, = rotations about r, 0, and z-axes. 

A dilation. 
A, p Lamé’s constants. 
Co stresses. 
a thickness of the slab. 
p = density. 
C z/a, the non-dimensional axial distance. 
¢, &, G5, C1, Cy = velocities. 

v, 0, Vv}, Vv, = non-dimensional velocities. 

p frequency xX 27. 


wavelength. 
y wave number. 
Vv Poisson’s ratio. 
Y a function of Poisson’s ratio. 


1, B = arbitrary constants. 


t time. 
mf Bessel functions of order n. 

T(r.2 t)) ‘ . 

; functions. 

g(r, z,t)} 
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/ 


h = a parameter = /( < fl —y?). 
N \A 


+2 


k = a parameter = [(ee—y') 
N 


bh 
8s ya. 
T pa. 
u v/s ha ya. 
v = y/s = ka/ya. 


1. Introduction 

THIS paper presents an elastic analysis of stress-wave propagation in an 
infinite slab of arbitrary thickness for the radially-symmetric modes under 
suitable line-source excitations. It then presents appropriate numerical 
data and discusses the results in terms of some recent experimental data, 
The treatment of cylindrical stress waves was first applied to rods by 
Pochhammer (1). Among later workers on this problem have been 
Love (2), Bancroft (3), and Davies (4, 5). The corresponding problem of 
stress waves in a plate shows more recent activity, notably by Goodman 
(6), Kane and Mindlin (7), and Feder (8). Correlatton between theory and 
experiment is difficult because of the complexity introduced into the prob- 
lem by the numerous boundaries of actual specimens and uncertainties 
regarding the exact form of the excitation. Also experiments on scabbing 
which necessitate the use of explosives to initiate stress waves, involve 
important non-linearities in the modes of propagation. This, however, does 
not invalidate the need for extending the ‘exact’ linear analysis, nor does 
it even reduce its importance, since it is necessary to ascertain the limits 
of its applicability. 

Goodman, in (6), has worked out the boundary-value problem for infinite 
plates and obtained the ‘frequency equation’, that is, the transcendental 
relation between frequency of vibration and wavelength. Kane and 
Mindlin’s paper gives an approximate solution for the high-frequency 
extensional vibrations of plates by use of energy methods. The still more 
general case of arbitrary line-source excitation can be dealt with by Fourier 
integral synthesis techniques based on the results of this paper, which should 
be of aid in experimental investigation of pulse propagation caused by 
explosives, such as were made by Feder and her group. 

2. The boundary-value problem 

Fig. 1 shows an infinite slab of thickness 2a excited along a line coinciding 
with the axis of z together with the appropriate cylindrical coordinate 
system. Details regarding the excitation will be discussed later. In this 
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shall restrict ourselves to a radially-symmetric 


squsoidal excitation. 











Fic. 1 
The displacements w,., wp, u satisfy the Pochhammer equations in the 


m given by Love (2), as follows: 


2 » 

¢ u \ a , ¢ A =f CW ») CWap ( l ) 
p (/ ZL) =p 

ct? C} r oo Oz 

( i | oA CW, CW. 

») » ? » = ) 

p A aL) =p =p (2) 

ot- roo Oz a] 

C-Uu m= oA 2u c 2u CW, 9 
p : A+ 2u) (7wa) ’ (0) 

ct- Cz Yr or r ob 


which A and » are Lamé’s constants for the material, p is the density, 


iA is the dilation, given by 


l c(ru.) Ll ou cou 
A ee (4) 
cr r oo Oz 
ithe con ponents ot rotation W,., Wa, W by 
| Ll ou cu ~ 
a. (5) 
2\r o6 Cz 
L/ou ou ' 
P | |; (6) 
2\ oz or 
l /e(rup) CU, = 
. | v ‘). (7) 
27 or 06 


\t this point the statement of the problem is mathematically similar to 


it lor the cylindrical rod—indeed, both are constructed from the same 


sic solution. For an angular frequency p and wavelength A, the only 


roduct-type solutions to Pochhammer’s equations are of the form 


U, A .e'42F (asr), u 


} 


A,e' 2 J,(ar), (3) 
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with x3 + ue is xa -5 (9o 
A+ 2p 
or tag = PP (9b 
be 


For the rod problem we obtain the required two linearly-independen; 


solutions by putting 
—s 
yy = y = 2n/A, 
v2 = h,k according as (9a) or (9b) holds. 
For the present problem of the slab, on the other hand, where thy 
boundaries go with the coordinate z instead of r, we obtain an appropriat 


pair of linearly-independent solutions by interchanging the roles of , ani 


¥, le. putting a, = y = 27/A, and in one solution 


92 9 
eo xe *) ai 


while in the other x4 k |= 1), (101 
AY be 


With «, and a, understood as above, and keeping in mind that (8) is 


¥ 


solution, the general solution may be written as 


u, = (A coshz+ Beos kz)J,(yr)e™, (1l 
U. - (-* sin hz OY ian kz) Jay -. (12 
: r i 


which can be verified by direct substitution in (1), (2), and (3). 
The boundary conditions that the front and back surfaces z = a ai 


z —a be free of stress are 


CU. CU,, u U. 
o,, = AA+2n— A—t+A—+(A-4 2u)< 7=—0, (12a 
a Oz or a 02 : 
‘Cu , 1 eu, 
0.9 H| — + : 0, (12b 
: oz r 00 
Cu, cu, 
cd. pi — + — 0. (12¢ 
, OZ cr 
When these are worked out and the common radial factors cancelled 
we obtain, for r + 0, 
ae 
A cos ha|Ay+(A+2u)—)— Beos ka 2py = 0, (I: 
vy 
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A sin ha(—2h)+ Bsin kal - ip 41. _ = 0. 
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these equations to be satisfied by non-identically-zero values for A and 
requires that the determinant formed by the coefficients of A and B 
zero : 


h? y , 
sha sin kaj Ay+ (A+ 2p) \ k4 t) +-sin ha cos ka(—2h)(2wy) = 0. (15) 


/ 


This relationship between p and y, or alternatively between the propa- 


ation velocity and wavelength, is referred to as the frequency equation. 


nother form of the above equation was first obtained by Goodman (6). 


1 Thin slabs 


As the thickness a of the slab tends to zero we can determine directly 


he limiting wave velocity of propagation. Given that the frequency p is 


nstant, two distinct states may occur: (i) either the velocity ¢ = p/y 
comes infinite, corresponding to y — 0, or (ii) the velocity is finite for 
y. The first corresponds to the second and higher modes and the latter 
responds to the lowest, i.e. the first mode of propagation. This mode 
nonly be propagated in the limit a > 0 provided y does not tend to 
fnity, so that # and k also remain finite. We then have 

cos ha ~~ |, cos ka ~ 1, 

sinha ~ ha, sinka ~ ka, 


ithe frequency equation (15) becomes 


h2 v2 ; 
ka Ay 1 (Xr 2) : \ k4 1 4ah* wy 0. (16) 
nce y = 0 we may divide by it, then solve the resulting equation for the 


tio py, which gives, for the limiting velocity, 


Pp * 9 (“) I( A TPH ). (17) 
Y N \p] N \A+2p 


r steel we have é 1-68c,, where c,  (u/p) is the velocity of distor- 


mal waves. This value for € can also be obtained by considering the 


topagation in a thin slab to be a plane stress problem. In plane strain 


e velocity of plane dilatational waves is , (A+ 2y/p), and so to obtain the 


responding velocity in a plane stress problem, replacing A by 2Au/(A+ 2) 


see Love (2), p. 208) we can get (17). 


1, Determination of the lowest mode 


By suitable transformation the frequency equation (15) can be put into 
orm which contains three variables only, thus facilitating the numerical 


mputation procedure: 


h2 v2 h tanha 
) S —4—p 
k’ tanka 
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This may be written further as 


2 2) ' 4h k tan ha 0. (1s 
yy tanka 


» 
C3 
where c¢ = wave velocity p/y, 
C, = distortional wave velocity ,/(j/p). 
: h I (c? \ k (ie* 
We also have iI . 1), = Jl = se 1), (19 
Y WN \G Y C3 
where c, is the dilatational wave velocity given by 
‘ A+2u 2u(1—v c 
eG = f I ) = fee (20 


p (1—2v)p ¥ 


and x2 


h / 2 : 
/ | =~ a*—1]}. (2 
Y N \C2 


c/e, = V, (22 
we have 
(V2—2)?4 4,/(V2x2—1),/(V2—1) aT - ih. (23 
When L< V2 < I/e*, 


equation (23) may be written as 
mts a oe = tanh ay,/(1— V2a2 
(V?—2)?—4,/(1— V2a?),/(V2—1) ayy | -s ST wi (24 
tan ay,/(V?—1) 


This equation is to be used for the calculation of the first mode. Nov, 


from (77), P 9 (") i A 4 E | — 2 
J p A A+ 2p : 1—v]’ 


so that V2 V2s2 (1—v)! (25 


fo - 
where V c/é. (26 


By giving various values to Poisson’s ratio v and solving for ay from (24) by 
giving various values to V, we locate the points on the curve of the first mode 
as given in Table 1 and shown in Fig. 2. In Table 1 the values have been 
computed so that the maximum error in the third decimal place does not 
exceed +3 units. For Poisson’s ratio v = 0 the point where V becomes 
unity is ay = 1-5708. Since the curve for v = 0 is to be understood only 
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in the limiting sense v > 0, we have chosen that branch which remains 


identically equal to unity for all ay less than the above value. 


TABLE | 


Values of V as a function of ay and Poisson’s ratio. for the lowest mode 


[’ 
I 2 3 1 o's 
i ) I I I oO 
197 93 pol 0°90 
7 54 )55 Ir | o-854 
12 355 793 | 0723 
SI 777 734 53 | 0°626 
714 "OSS 054 C 14 500 
65 "6039 61 73 0°52 
71 2 | 0°592 548 | 0-500 
1-0 <= ne 
0-9 - t - + — 
v=05 0-4NO3NO0:2N01A\\0-0 
0-8 
r 
V ; 
( 
0-7 ts-4 
0-6 = 
“al 
0-5 => 
0 1-0 2:0 3-0 00 
S=ay =27a/A 


Fic. 2 

The short-wave limit, that is ay > «, can only be obtained from (23) or 
24) if J 
the numerator and denominator of (23) become hyperbolic and we have 
V 20) 


itanhay, (1—V?) > 2, 


> becomes less than unity. For then the tangent functions in both 


tan ay, (V2a?—1) = itanhay,/(1 >i, 
tan ay, (V*—1) 

ind the frequency equation reduces to 
(V?—2)*—4,/(1 


or, putting V? = y, 
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This equation is the same as that deduced by Rayleigh for the velocity of 
surface waves (Rayleigh waves) on a semi-infinite isotropic solid. This 
equation has always just one real positive root, which gives us the asymp- 
totic value for the curves, i.e. the surface velocity: 

C, = 2980 m/sec (steel). 


For the second and higher modes, (23) can be used, but the arguments 
of the tangents become very large so that the location of the roots becomes 
quite cumbersome to calculate. For the determination of these modes a 
method described in the next section has been used. 


4-0 














5. Determination of the higher modes 

The higher modes become of importance as the frequency of the excita- 
tion increases, or if the distribution of the exciting forces along the axis of 
the slab changes rapidly. The method described here was used to calculate 
the second and third modes explicitly, and can easily be extended to modes 
of higher order. 

If, in (15), we introduce the variables 

xz = ha, 8 = ya, 


y = ka, T = pa, 
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the frequency equation, now independent of the thickness a, becomes 


v2] / 82 x. 
}A+(A+ 2p) AI . | cos sin y 4. —sIn & COS Y Q, (28) 
| lly ; 
, pt ~ OF ; 
rith x ——s yy? — 2, (29) 
A op pL 


‘quations (28), (29) form a system of simultaneous equations for the 


wiables x, y, 8s, and + 


We seek ply | 7/8). 
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Fic. 4 


\ direct substitution of the relations (29) into (28) would lead to a 
rmidable equation to solve. Instead, if we eliminate 7? from (29), we 


tain (A 2)a2 py (A | p.)s® 0. 


ividing by s* and letting 
ris u. Y/s v, 
e obtain a relation which is fixed for a given material 
(A 2)u* pe (A+ jt). (30) 
is relation is shown in Fig. 3 for steel. The dotted branch corresponds 


imaginary values of wu. The frequency equation (28) now reduces to 


(v2— 1)? cosxasin y+4uv sina cosy = 0. 
he final changes of variable 
(ety) =& (@-y)=y 
luce this further to 3 / sin €/sin », (51) 
here \ Ii (v2—1)?+ 4uv], 
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4uv}. 
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Since v is a known function of uw, as given by (30), these are also functions 
of uw. In this analysis it is convenient to take wu as the independent variable, 
We finally have (31) and the relation 


x/y w/v q(u), (; 


rt 


which is now in a form suitable for fairly rapid graphical solution. 

Fig. 4 shows some of the loci 8/« = const, as given by equation (31), in 
the (€, 7) plane, but labelled by the corresponding value of wv. It also shows 
the straight lines w/v const, labelled in terms of the same values of 
The intersection of the two sets gives the common solution (uw, €, 7), fron 
which we go back and determine x, y, s, and 7, then p and y for a given ¢ 
and hence finally the velocity c. 


Steel 
& =5400 metres /sec. 


Cs=2980 metres /sec. | 


3-0 


>) 


2:0 








1-0 2:0 3-0 
s=ay =27ra/A 
Fic. 5 


The main numerical results are shown in Fig. 5, which plots the ratio o! 


was chosen to be similar to the curves which Bancroft (3) obtained for the 
cylindrical rod, to afford easy comparison. 

While the method of this section does not lend itself easily to the lowest 
mode, there is one point on that curve which is easily obtained. If we make 


v 1 then (30) also gives for steel u 0-6407 from which x 1-048), 
y = 4x, s = 4x, c = 4580 m/sec, or c/é = 0-85. This corresponds to the 


point marked A in Fig. 5 which coincides with the value obtained from 
formula (24). 


ce to €in terms of s = ay for the first three modes. This manner of plotting 
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6. The excitation 


In this section will be discussed the conditions under which a wave repre- 
sented by the solution (11), (12) can be generated. Since the boundaries 
have been made free, there can be no excitation from the faces of the slab. 
The only remaining source for the wave, which must be radially symmetrical, 
can be a line corresponding to the axis of the cylindrical coordinate system 
or on a cylinder concentric with this axis. 

We have, from (11), (12), (12 ¢), 


[—2hA sin hz+ B(—k+y?/k)sin kz] J,(yr)e™ (33) 


nd 0, AA 2p 


J (yr 
(A cos hz Beos iz) (A Qu)‘ (yr) ! A 
cr 


—J,(yr) eipt (34) 
r 


O.0 0 


for the values of the stress components on a cylinder of radius r. In these 
expressions the constants A and B are no longer both arbitrary, but must 
be connected by the boundary relations (13) or, equivalently, (14). 
We have, from (13), 
A 2uy cos ka 
B COs ha\Ay (A+ 2u)h* y] 


cos ka 2uy coska 2 (35) 
35 
cos ha py(c?—2c3) cos ha V2 


The case of axial excitation can be dealt with by letting r— 0 in these 
expressions. 
From (33) o,. (} so that only the radial normal stress of the above is 
present in this case. Its value then is 
Cleoskacoshz 2 


- cos kz | e', (36) 
cos ha y2—2 





where C is an arbitrary constant replacing 26(A+-2u+-1). In terms of new 


variables this can be re-written as 





, Clcosy 2 tae x oo 
(F+r)r=0 p(C) - “= cos x¢ + cos yf je” a) 
y |cosx V?—2 
(C/y) Fe, say, 
where { = z/a varies from 0 to +1 as z varies to +a. Fig. 6 shows how 


the normal stress varies along the axis in the case of steel, for some differing 
values of the parameter s or 7, for the lowest mode, J. 
From a physical point of view, that curve of the figure (and with it the 
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corresponding value of the parameters) has to be chosen which most nearly 
corresponds to the excitation. If we imagine this to be caused by an explosive 
packed into a very thin hole, then the curves of normal stresses shown in 

















the figure, together with the fact that o,, = 0, fits in with the picture of 

as | i i a. 

| s=0 
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= | a 2 ee 
| 7-083 
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an expanding gas at high pressure. It must be borne in mind, however, 
that our present analysis only applies strictly to a periodic disturbance. 
Because of the Fourier Integral Theorem, it is possible to extend our 
inference to pulse propagation. 

Although little is actually known at present concerning how the above 
pressure varies along the hole, we can reasonably imagine that due to 
some pressure relief it will be expected to fall off towards the ends some- 
what. The one limiting case of complete confinement (pressure constant) 
is that for s = 0. The other is that for complete release at the ends, corre- 
sponding to the value s = 1-4 (approximately). It is interesting to note 
that beyond this value there is a change in sign of the stress to tension 
over part of the range. This, of course, could not be excited by means of 
expanding gases. 

The conditions of continuity of the medium in the case of a radially 
symmetrical disturbance necessitate that the radial displacements at the 
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CYLINDRI( 





line of symmetry must be zero. This is evidently satisfied in (12) above, 


since U, Oatr 0, 


7. Experimental results and their comparison with theory 

Having completed the general theoretical analysis of the infinite plate 
we wish to discuss here from its standpoint two recently obtained experi- 
mental data. Since the data must of necessity be limited to finite plates, 
it can be expected to be an approximation to the results of the theory for 
infinite plates at best. However, it affords some indication which can serve 
to illuminate both approaches. 

Some work was done at Bell Telephone Laboratories on the propagation 
of longitudinal stress waves in circular disks of barium titanate with various 
diameters, thicknesses, and their ratios, i.e. d/2h R/h, where d is diameter 
ind R is the radius of the disk. The resonant frequencies of the disk for 
the first two modes were determined with the help of a tuneable oscillator 
and a frequency counter. The results of the above work were first used by 
Kane and Mindlin (7) and the same have been used here for calculations 
in Table 2. 

As also indicated by Kane and Mindlin the frequencies of the first mode 
of vibration were comparatively easy to determine but the measurements 
forthe second mode were difficult because in addition to locating resonances, 
their relative strengths had also to be determined. Thus the frequencies 
of the second mode are less reliable than the first. 

In the simplest way that resonant radial vibrations of a disk may be 
excited, the disk has a node at the centre and the next antinode at the free 
periphery. This is also the easiest vibration to excite in experiments. Thus, 


taking indiach (37) 
the velocity of propagation c becomes 
Cc Af da /2z7. (38) 

Also ay 27(a/A) 2n(a/2R) = x/(R/a). (39) 
The elastic constants of the barium titanate used for the disk are 

1 0-3 

He 49-3 x 101° dynes/cm?, 

p = 5°33 g/em'., 


Hence the thin plate velocity is 
é y (MH) p)/{2/(1—v)} = 3635 m/sec. (40) 


With the help of (37), (38), (39), and (40) values of c/é and ay were calculated 
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and are given in Table 2. The dimensions of the disk and corresponding 
resonant frequencies are taken directly from Table 2 of Kane and Mindlin’s 
paper (7). om 

TABLE 2 


Velocities calculated from experimental data for barium titanate disks 











2 
(v 0-3) 
Diamete Thickness 
(cm) (cm) Ist mode 2nd mode r 
—__—_—— oes, i 
2k 2a cli cli R/a 
2°05 0°045 0°5027 2°2225 o°755 
2 j / 
2°63 0°640 0°8927 2°2176 07644 
2°16 0°564 orsg2%5 1007 O'OIO! 
2°16 0°037 Osgig 2°0525 0°9207 
2°17 0°638 0°8935 2°0044 0°9207 
1°58 0°596 o"gos4 1°9037 I°1S1o 
I°s9 o'612 o'gl2 1°5459 I* 2003 
rs8 604 O°907 1°SQ77 L°19Q9I 
I°31 0°639 05504 1°4795 1°5325 
1°OS 636 O°5245 I°2421 I*go4 
1°05 0°637 #5349 1°2330 I°Q15 
783 0°629 o'7 194 1°2322 | 2°5122 
| 

ree of => as ‘ete 
OTdi 630 7230 1°2377 2°5751 


The experimental points of Table 2 are shown in Fig. 5. The theoretical 


curves drawn there are for steel (v 
0-3). 


to the different Poisson’s ratios. Bearing in mind the fact that the experi- 


()-29) but the experimental points are 
for barium titanate (v Part of the discrepancy in the two is due 
mental points were obtained for finite disks and the curves for the infinite 
plates, the agreement between the two is fairly good, especially for the first 
mode. 

It is worth while to discuss here two bases used for doing the above 
calculations: 

(i) The velocity of propagation of the cylindrical stress waves in infinite 
flat slabs can be determined from experiments with circular disks of 
finite radius, only as an approximation. — 

(ii) It is assumed in the above calculations of both the first and second 
modes, that under resonance the centre of the disk is a node and the 
periphery an antinode, there being no node or antinode in between (cf. (37)). 
This is evidently true for the first mode but is also true for the second or 
higher modes. Otherwise the velocities for the second mode would be lower 
than the velocities for the first mode, which does not make much sense 
and hence the above assumption can be justified. Actually, the first mode 
differs from the second in the number of nodal planes of the lateral displace- 
ments u,. The first mode has one nodal plane, the second has two, and so on. 
In a second group of experiments made at Frankford Arsenal by Feder 
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tal. (8), a stress wave was initiated in a transparent plastic (Columbia 
Resin CR—39) by dynamic loading methods. In order to ensure a wavelength 


sufficiently short compared with the specimen and for other reasons a 


lasting cap was exploded at a point on the specimen and photographed 
photoelastically with a high-speed framing camera. Among the variously- 


shaped specimens one of them was a circular model as shown in Fig. 7. 


A comparison of the relevant velocities is given in the following table: 


1. Measured fringe velocity (in./microsec) 0-104 

2. Caleulated velocity of dilatation waves (compression) 0-077 

3. Long thin bar velocity 0-0414 

1, Thin plate velocity 0-051 

5. ‘Thick’ plate velocity 0-027—-0-051 


In the reference it is noted that the measured velocity is greater than the 
ngitudinal velocity c,, by a factor of 25 per cent. The non-linearity of the 
stress-strain relation at high stress would not be able to account for this. 
\ctually it would be more accurate to compare the measured velocity with 
the ‘thick-plate’ value, i.e. that given by the analysis of this paper. It can 
seen that this would furnish at most the upper limit € = 0-051 in./micro- 
vec. The pulse length, although not accurately known, is sufficiently short 
tomake the assumption of infinite plates reasonably accurate. This can be 
letermined from the photographs which were taken in connexion with these 
experiments. It is not possible at the present moment to determine the 
elocity more exactly because of this uncertainty, and because, as pointed 
ut by Davies (5), the detonation phenomena at the interface between a 
solid and a gas are too complex. 


Nevertheless, the higher observed velocities depart even more by com- 


| paring with the plate velocity instead of the dilatation velocity. It may 


be due to the presence of shock waves generated by the high degree of 
ading used. Another possibility may be that second or higher modes are 


eing obser V ed. 
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8. Group velocities for the first mode of propagation 
The relation between the group velocity c, and the phase velocity c is 
A (de/dA), 


which on dividing by the thin plate velocity ¢ can be written as 


Cy c 


A V—A(dV /dA), 
or, putting A in terms of s, 
A V +s (dV /ds). (41 


TABLE 3 


Group velocities V, against ay for various Poisson’s ratio 





Ve 

v 
ay o'o orl 2 03 o'4 O°5 
fo) 1000 1*000 I 1°00 1°OO 1*000 
o°5 1°000 | 0°999 “O91 0°973 | 0°939 | 0°883 
O I‘000 “986 32 | 0°839 | 0°723 | O°592 
cS 1*°000 0°734 0°568 0°477 | 406 0°375 
0°233 | 0°434 $71 23 373 | 07352 
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This formula is used for determining the group velocity curves of : 
gainst s, graphically, from the phase velocity curves of V against s. The 
phase velocity curves were plotted accurately on a large scale and the values 
ftangents read. This gave dV /ds and consequently V, against s. The error 
ntroduced due to inaccurate reading of the tangents was reduced by 
lotting curves of V, against Poisson’s ratio, v, for various ay and smooth- 
ng them. These curves were quite flat and so minor errors in points could 
ily be detected. The results thus obtained are given in Table 3 for 
ious Poisson’s ratios and plotted in Fig. 8. Due to the graphical proce- 
lure, the third decimal place in Table 3 may not be accurate. Itis interesting 
note that all the curves appear to pass through ay = 2-2, V, = 0-325. 
\s in the case of rods, the group velocity V, approaches the surface velocity 


r very small wavelength. 
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SOME CONDUCTION PROBLEMS IN THE HEATING 
OF SMALL AREAS ON LARGE SOLIDS 


By P. H. THOMAS (D.S.I.R. and Fire Offices’ Committee, 
5 Fire Research Station) 


[Received 19 September 1957] 


SUMMARY 
This paper discusses a number of problems the common feature of which is th 
heating by a constant flux of small areas, generally circular, on the surface of a 
semi-infinite solid or on an infinite thin sheet. In addition to a number of analyti¢ 
results, the variation with time of the surface temperature of the centre of the heated 
area on a semi-infinite solid has been computed for various sizes of area and degrees 


of cooling. 


1. Notation 


x,y,z cartesian coordinates. 


r,« cylindrical coordinates. 

iy temperature above ambient. 

t time of heating. 

K thermal conductivity. 

p density. 

c specific heat. 

k thermal diffusivity (K /pc). 

Q incident heat flux per unit area. 

R radius of heated circular area. 

2L width of heated infinite strip. 

D perpendicular distance between point source and surface of solid. 
A thickness of thin sheet. 

H heat transfer coefficient at the surface. 
h H/K. 

V2 Laplacian operator. 


2. Introduction 

WHEN a small area on the surface of a large solid is heated, the temperature 
at first rises in substantially the same manner as does that of a much larger 
area. Eventually, however, the cooling effect of heat flowing laterally away 
from the heated area becomes important and it is necessary in experimental 
work which employs small areas to simulate larger ones to be able to assess 
the time before this occurs. For example it is desirable to know the time 
after which the temperature of the centre of the heated area is no longet 
within, say, 10 or 20 per cent of what it would have been had the heated 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 4 (1957)] 


wea EX 
ne-dil 
the me 
Con 
solid, : 
oss fr 
the rat 
two lo 
oss pl 
repres 
ver t 
f hea 
numb 
neglig 
may | 
kind 1 
if its 
Iti 
f sm 
given 
creul 
when 
temp 
i) als 
rea | 
is QV 
const 
ture 
sh 
Ty 
of ras 
over 
of a 
prob 


equa 


The 


and 





NG 


olid. 


menta 


» ASSES 


ie time 


longel 


heat 














SOME CONDUCTION PROBLEMS 483 


rea extended over the whole surface, i.e. had been a plane source producing 
ne-dimensional flow. An approximate assessment can often be made from 
he magnitude of certain dimensionless parameters. 

Consider a small circular area of radius R on the surface of a semi-infinite 
olid, and suppose that only this area is directly heated. The rate of heat 
ss from unit area of surface may be taken to be of the form H@ whilst 
he rate of heat loss by conduction into the solid is K@/R. The ratio of these 
wo losses is H R/K and if this is very much larger than unity, the surface 
iss predominates and the size of the area is unimportant. Similarly kt/R? 
epresents the ratio of the rate of lateral loss of heat A@/R by conduction 
ver the cylindrical surface of a region of unit depth to the rate of retention 
fheat Rocé/t by that region. It follows that if the ratio kt/R? (the Fourier 
number) is markedly different from unity then one or other factor is of 
egligible consequence and an approximate but simpler conduction model 
lay be substituted. Where kt/R? ~ 1 and HR/K ~ 1 analysis of some 
‘ind is required, if it is desired to use an experimental system to the limit 
f its accuracy in representing large areas. 

It is the purpose of this paper to give a number of results for the heating 
fsmall parts of large solids or thin sheets. In section 3 some results are 
siven for the transient and steady-state temperature at the centre of a 
ircular area on the surface of a semi-infinite solid subject to surface cooling 
when the area is heated by uniform constant heat flux. The steady-state 
temperature on the surface of a large solid due to an external point source 
salso discussed. Section 4 deals with the problem of a heated circular 

rea or infinite strip on a thin sheet. Finally, in section 5, a brief discussion 
sgiven on the results obtained in sections 3 and 4. Note that all physical 
onstants of the heated material are taken to be independent of tempera- 
ture variation. 

3, Local heating on semi-infinite solid 

Two conditions are considered. The first is the heating of a circular area 
fradius R on the surface of a semi-infinite solid by a constant flux uniform 


ver that area and zero elsewhere. The second is the heating of the surface 


fa semi-infinite solid by a point source external to it. In this latter 
problem only the steady state is considered. The transient heat conduction 


equat ion is 


ro = &. (1) 
c 
The initial condition is iy Q, t 0 
ind one boundary condition is 
6 —> 0. x —> OO. 
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When the heat flux is incident on a circular area, the second boundary 
condition is that at « = 0 
06 “.* é - 0 2; 
He—K*® = Y (r<R and ¢t ); (2a) 
ex |\0 (r>R and t>0). (2b) 
For the point source external to the solid it is 
- 00 
Ho—K— = Q(r) (t> 0), 


C2 

where Q(r) can be determined as a function of r. 

In practical problems, the primary interest is in the temperature of the 
centre of the heated area. 

(i) Transient temperature at the centre for constant heat flux 

The transient temperature at the centre of the heated area may be found 
by making use of Green’s function u for a point source of unit strength at 
(x’, y’,z’)in a semi-infinite solid having surface loss of heat at the boundary 


«x = 0 (see Carslaw and Jaeger (1)). The point source is considered to be 
situated at the surface x’ = 0 so that the transient temperature at the 
centre of the heated area (x1 = y = z = 0) is 
Ug = ———, e424 _ . erfe hy (kt)exp me 2 - (3) 
4(akt)? tk ’ 4kt 
Hence for a continuous source of strength Q/pc, over the region r < R, 
where r , (2’*+-y’*) the temperature is 
t R 
Qf "fo nee _ 
0 | dx | ~— eh*kA erfc h,/ (kd) \e-7*/4*r dr. 4) 
a 7 oe J \at(kayi ka v{ \ 
0 0 


Integrating first with respect to 7, writing 0, for the temperature with 
no surface cooling and substituting A = w?/h?k we have from equation (4) 





hv (kt) 
20 ake cade i 
4, Ove— H | (l—e h® R*/4@ \we” erfe w dw (5) 
0 
QR 2 (kt\4 . R 
and Ove : 2(] —e—R*4kt) 1 erfe mh 
un NC K ar (l—< )+erfc =| 


a result obtained by Blok (2). 

The integral in equation (5) has been numerically evaluated and the 
resulting values of 6, are shown in Fig. 1, for various values of AR and 
kt/R*. It is possible to find the values of kt/R? at which the temperature 
at the centre of the disk is, say, 10 or 20 per cent less than it would be for 
a source extending over the whole surface x. = 0, i.e. a plane source. These 
values of kt/R* are given in Fig. 2 as a function of hR, while in Fig. 3 the 
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corresponding values of h?kt are given as a function of hk. For any h, R, 
and k it is thus possible to find the maximum duration of heating for which 





the temperature of a small area is sensibly the same as a large one when 
heated by a constant uniform flux. 
20% difference 10% difference 


h2kt Asymptotic tohR=4-6 Asymptotic to AR- 0 


| 
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Fig. 3. Values of h2kt as a function of hR for 10% and 

20% difference between centre of circular area and 





extended plane 


(ii) Steady-state temperature at centre for constant heat flux 
For this problem a solution in terms of tabulated functions is possible. 
A solution to equations (1) and (2) for the steady state is 


QR [ e-%*J,(ar)J,(aR) 
a h-+-a 


0 


6 dx, (6 


where J, and J, are the Bessel functions of the first kind of zero and first 
order respectively. 

This integral can be evaluated at r = 0 and x = 0 by means of an 
integral given by Watson (3a) to give 


QR ; l 


K "AR b7{¥(hAR)—H,(AR)} |, 4) 
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shere Y, is Weber’s Bessel function of the second kind of first order and 
f, is Struve’s function of first order. The variation of K6,/QR with AR is 


hjown in Fig. 4. As hR >, 6, > (QR/K)(1/AR), ie. Q/H, which is the 


0 


teady-state temperature for an infinite plane source—as should be the 


ise. Furthermore, when h is zero, i.e. there is no surface cooling, it can 
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Fic. 4. The steady temperature of the centre of a heated circular area 


on a semi-infinite solid or below a point source 


shown from equation (7) that 4, is equal to QR/K. This case, however, 
sof special interest since 6 for all values of r may be obtained in terms of 
tabulated functions. 

By considering sources of strength 7, dr,d¢, Q/pc over the region 


ry; R 0< d < 2zn, x 0, 


the steady-state temperature distribution at any point on the surface x = 0 
in be obtained as 


Qf, - i; 
A "ee | dx | dd | rexp| a 


0 0 0 


2 | 42 Ire. CCAS dA) 
"7 i) 2rr, cosh (kd) ? dry. (8) 
4ky 
We make use of a particular case of Weber’s second exponential integral 
3b), namely, 


) ap x 
. 


l | exp| - "I a3 sone dd = | e-vka F ( xr) Jol xr Jou da. 
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After integrating with respect to A and r, we obtain 


= 5 , (of) Jo(ar) dx. 
; mJ x 
0 
It may be shown (3c) that this reduces to 
9 
9 — 2 OF rR) (r< R), (9) 
a K 


-2$le()-(-S}eG]] mm 
a K r uz ij 


where K and E£ are the complete elliptic integrals of the first and second 
kinds respectively. 
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Fic. 5. Surface temperature distribution due to uniform circular area over 
0 r < FR and to an external point source 
This distribution of temperatures is shown in Fig. 5. In particular the 
maximum temperature at the centre (r= 0) is QR/K and the mean 
temperature over the area r < RF is 8QR/3rK. 
(iii) Steady-state temperature due to an external point source 
O 


™ 


P. 
| 
DD | 





Fic. 6. Point source external to semi-infinite solid 


If W is the heat radiated per second at O a distance D from a non- 
reflecting semi-infinite solid, the heat flux at P normal to the surface of 
the solid is (see Fig. 6) 

W cos 6 WD 


Oe’ on some nema 
wl") = F(D?-pr8) ~ Fx(D=p7)! 
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. dé WD 


herefore, at x 0, hé Se: lla 
1erefo z a inK(D®47)! ( ) 
nd at x 0: 6 0. (11 b) 
\pplying the Hankel transformation 
8(p) [ A(r)rJ,( pr) dr, (12a) 
0 
6(r) ( O(p)pJy(pr) dp (12 b) 
0 
ve have for the steady-state from equations (11a), (11 b), and (12a) that 
}— Ae-P? and at x = 0 = 
(ht pb =o [ ler) ee _ TT ee (13) 
torK J (D?+r?)! 4K 


0 


from an integral given by Watson (3d). At the centrer = 0, from equations 
12b) and (13), 


O a 1—hDe"P | —— : (14) 
4nK D : u 
hD 


The integral in equation (14) is the exponential integral and is a tabulated 
function. The value of 6, for various values of hD is shown in Fig. 4, 
) being denoted by Rk and W/47D? by Q for comparison with equation (7). 
For no cooling, equations (12 b) and (13) give 


we W 
@ B ‘\e—PD d apomeeneteaaiat 
ink | “\Rrie’ “P = CK (D*+P) 


(15) 


sing an integral given by Watson (3e). 

The central temperature without cooling is seen to equal that due to a 
lisk on the surface of radius D and intensity equal to the maximum 
V/47D?. Because the total heat falling on the surface x = 0 is 27D? 
times the maximum intensity, and for a uniform disk on the surface it is 
nly 7D? times the maximum, there is a factor of 2 between the respective 
temperatures at large radial distances. 


4. The heating of thin sheets 
In these problems the thickness A of the thin sheet is assumed small 
enough for thermal gradients across the sheet to be neglected. The validity 


f this assumption depends on the conditions 


hA A<R, or AKL. 





l, 
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In transient problems the equations will only be valid for 
t > A?/k. 


The systems considered are the heating by a constant flux of a small 
circular area and an infinite strip as shown in Fig. 7. 
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Fic. 7. Heating of thin sheet 





(i) Circular area 
For the condition shown in Fig. 7 we have 


220 100 120. 2h0 Q(r) 


| | : 16 
or? rér kat A KA 
( ad 
where Q(r) wd : me (16a 


6=0 att =— 0. 
Applying the Hankel transformation to equations (16) and (16a) we obtain 


29 1d6 QRJ(pR) , 2hé 
)* ~-—— “| zs 
: ik di KAp A 
whence ; 
QR [ J(pR)J(pr) , . 
-| 1—exp{ 21. 2h/A)kt}] dp. 17 
KK } “(pp ahjdy 11 exPt— (Pt 2A] dp 
0 
For no cooling the transient central temperature is, from equations 
(12 b) and (17), 


L 


6, Qt ( o-X*4 X2 | e-u “), (18) 
pcA u 


. 


Xe 

; R 

where X : : 
2, (kt) 

The relation between 6, and X is shown in Fig. 8. The steady-state central 

temperature with surface cooling is, from equation (17), 


OR : Be a dp, 
KA } (p?+2h/A) 
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6, 
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to be 


which, from an integral given by Watson (3d), may be shown 


Q | 2h\ .. e } 
Qf) rp VIR 
"0 = o7\ JG) Ja 


(19) 








pcg, 
Qt : 
] Str 
6 
4 | 
‘ | 
| 
0 02 04 06 O8 10 12 14 #%16 18 20 
R — 
2V (kt) 2V(kt) 
Fic. 8. The effect of size and heating time on the maximum tempera- 
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where K, is the modified Bessel function of the second kind of first order. 


[he relation between 6, and 2h R2/A is shown in Fig. 9. 
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(ii) Infinite strip 
We have 


oe 1 00 2h0 _ Q(y) (20a) 
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Applying the cosine transformation 


Ap) = | Oy) cos py dy, (21a) 
0 
ot 

Hy) = - | O(p)cos py dp (21b) 
7 


. 


0 
we obtain from equations (20a), (20b), and (21a) 
1d0 2h0 QsinpL 
kdt' A KAp ’ 





—p§ — 
and eventually 


sin pL 


9 ‘ 
6 > = - @ o | 
a KA J p(p?+2h/A) 
0 








[1—exp{—(p?+2h/A)kt}|cos py dp. (22) 


If there is no surface cooling the temperature along the axis of the strip is, 
from equation (22), 


( - 
ee (ext 2Y*erfe Y + — Ye-} }, (23) 
pcA N77 
: L 
where . 
2,/ (kt) 


The variation of 6, with Y is shown in Fig. 8. The transient axial tem- 
perature with cooling is given by 
t 


Q [ DQhkA/A f L ) 
- _v_ | ¢-2hk\A ey 1A, 24 
ped , ; ™ 2. (kA) ‘ (*4) 


. 


) 

so the effect of cooling may be estimated from the magnitude of 2Ht/pcA. 
The steady-state temperature with cooling is obtained from equation 

(22) with t> o and y = 0. Thus 


A . 
= (2 _ p—Lv(2h 4)), (2: 


bo 


which is shown in Fig. 9. 


5. Discussion 

A number of results are listed in Table 1. If, when comparing the disk 
source on a thin sheet and on a semi-infinite solid, the semi-infinite 
solid is regarded as approximating to a ‘thin’ sheet of thickness 2R, 
the conditions for a 10 per cent difference between a small and large area 
are the same. By analogy with the transient behaviour of a disk on a 
semi-infinite solid we should expect the permitted value of kt/R? (for a 10 
per cent difference) to increase with )R?/A and the shape of the curve 
relating hR?/A to kt/R* (for a 10 per cent difference) to be similar to that 
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given in Fig. 2, hR being replaced by 2h R?/A. The minimum value of kt/ R? 
for a thin sheet is, however, 25 per cent less (0-21 instead of 0-28) than for 


, semi-infinite solid. 


TABLE 1. Values of size parameters at which size has less than 


10 and 20 per cent effect 











Condition for less than 
10 per cent diffe rence 


Condition for less than 
20 per cent difference 








Problen Condit from plane at centre from plane at centre 
int source near | Steady-state with hD 17 hAD 7-5 
emi-infinite solid surface cooling, 
uniform = steady 
flux 
Disk on semi-in- | Steady-state with hR 10 hAR> 45 
finite solid surface cooling, 
uniiorm steady 
flux 
Disk on thin sheet | Steady-state with hR2/A 0 hR2/A 2-9 
surface cooling, 
uniform steady 
flux 
Strip on thin sheet | Steady-state with hL?/A -6 hAL?/A > 1-4 
surface cooling, 
uniform steady | 
flux 
Disk on semi-in- | Transient with | AR 0 kt/R? < 0-28 | kR 0, kt/R? < 0-52 
finite solid surface cooling, 2 0-35 2, kt/R* < 1-4 
uniform steady 4 0-65 
flux 
Disk on thin sheet | T’ransient, no cool- kt/R? 0-21 kt/R? < 0-32 
ing 
Strip on thin sheet Transient, no cool- kt/L? - 0-36 kt/L? - 0-65 


ing 
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ON FREE MOTION IN THE GRAVITATIONAL FIELD 
OF THE EARTH 


By G. C. SCORGIE (Calcot, Nr. Reading) 
| Received 27 September 1956} 


SUMMARY 


A literal approximate solution is given for the motion of a particle in the earth’s 
gravitational field as modified by ellipticity of the figure of the earth. In particular, 
the resulting torsion of the path is calculated. 


1. Introduction 

In connexion with the launching of artificial earth satellites during the 
International Geophysical Year, and of intercontinental ballistic missiles, 
it is of interest to investigate the motion of a particle in the earth’s gravita- 
tional field. In both cases a body is to be rapidly propelled into an orbit 
which is then traversed under the action of gravity and disturbing in- 
fluences of various kinds, air resistance accounting, presumably, for most 
of the disturbance. The present note is concerned with the solution of the 
basic problem of the free motion of a particle in the earth’s gravitational 
field when account is. taken of the ellipticity of the figure of the earth, that 
is, of the departure of the external field of gravity from a force varying 
inversely as the square of the distance from a fixed centre. The departure 
being small, the most suitable method of dealing with it is likely to be by 
a numerical perturbation of the ordinary Kepler motion. However, for the 
discussion of general characteristics of the motion, some interest attaches 
to an approximate literal treatment described here. 


2. Approximate solution 

3eing interested in the purely gravitational field, excluding centrifugal 
effects caused by the rotation of the earth, we take a non-rotating reference 
frame with origin at the centre of the earth, and choose spherical polar 
coordinates q" (n 1, 2, 3); q' is the radius, q? is the angle measured fron 
the earth’s polar axis drawn towards the north pole, and gq? is the angle 
measured about that axis. Since, in discussing the geometry of the path, 
we wish to use tensor notation, it is used throughout, with indices running 
from 1 to 3. 


{Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 4 (1957)] 
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The purely gravitational potential, the negative of the potential energy 


funit mass, at a point outside the earth is (1) 


| Yo+% 
; a 
vyhere V, 1 
q 
ab(\ 3 cos*q*) : 
Vo : , approximately, 
1\3 . 
(q ) 
a 3-986 x 107° em? see~?, 
nd b 2-225 « 10" em?. 


may be regarded as perturbing the Kepler motion corresponding to J). 
\lternatively we can write 
V V+, 
; ; a ab(1—3 cos?q") 
vnere V; 7 - 
q' R(q')* 
1—q!/ Ryab(l -3 cos*q") 
(qi)3 . 
nd R is a constant length comparable with the mean value of q! over an 
rbit. This procedure is useful because v, is a smaller perturbation than vp, 
nd the motion in the field }, can be found by quadratures. It follows that 


nto 1M | 1—2q1/3R| and —_ ore l—q/R\. 

qT od he whe os 
Thus, for example, pro. ided that qi does not differ from R by more than 
) per cent of R, the field J}, contains at’ least 95 per cent of the perturbing 
transverse acceleration caused by ellipticity of the earth’s figure, and at 
east 63 per cent of the perturbing radial acceleration. 

The problem of motion in the field V; can be reduced to quadratures by 
standard methods (2). Since g° is an ignorable coordinate the conjugate 
nomentum ps, is a constant m, say, and we have a system of the fourth 
rder, the Hamiltonian for unit mass being 


(1) 


2 ,{ m* 2ab(1—3 cos*q*)) a 
1" 


l > { ; 
H | (p,)2 Pe ; (q*) ” = 
| q 


2 (q')? \sin*q? Rk 

Since the Hamilton—Jacobi equation is separable, it is possible to find a 
contact transformation to new variables in which the coordinates and 
momenta are constants. Alternatively we can make an extended contact 
transformation, transforming the time while preserving invariant the 


Pfaffian 


A p,dq p,dq* H dt. (2) 


ind solve the canonical equations in the transformed variables. Adopting 
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the latter course, we regard the negative of the Hamiltonian as the momen. 
tum canonically conjugate to the time, noting that, since time does not 
appear explicitly in the Hamiltonian, its value must be a constant: and. 
since the coordinate system is inertial, this constant is the total energy £, 
say, of unit mass. Hence we have 


where H is given by (1). 
Choosing q' as the transformed time we can write 
A= —FEdt4 p,d¢* 


where the transformed Hamiltonian J is the value of - p, found from (3 
Thus 





Jdq', (4 


y= — [om (0 _ (qu afm 2ab(1—Beosty?)) 2a) 
(q?)* \sin?q? R a q 
From the canonical equations 
oJ _— dit eJ dq oJ dk od dps 
. CoE dq’ Cp, dq et = dq’ eq dq?’ 
we obtain ¢ and q? in terms of q', and q? in terms of q?, as follows. 
lq? . 
From . - 
dq’ J(¢')? 
ai dp, _ __ cos q°| hal ais 6ab sin q°| 
dq J (q¢ )?|sin3q? R | 
Ip, — 2 6ab ain at 
a owe dp, _cosq*{ m* _ 6absing?) 


dg —s py \sin'g?' = | 


which can be integrated to 








{ m*  6absin®g? \4 ; 
Pe — ees er er es ~ F T n =" (5 
| sin?q? R | 
where n is a constant. 
dt ee 
Hence =i q{2E(¢q)?+ 2aq!—(n-+4ab/ R)}-! (6) 
Cc q 
oi dq? _ {n—(m/sinq?)*+ 6(ab/ R)sin2q?}! (7) 
: dg @{2.E(q)?4 2aqi—(n +-4ab/R)\" . 
— : 1q? 
By definition, m = p, = (q'sin gy; 
C 
dq? m m? 6ab sin2q?) -1 
hence : s os a = ME cin nY va (8) 
dq? sin®q?| sin?q? R | 


The solution by quadratures is expressed by (6), (7), and (8). 
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t—t, = a(—2E)-?(s+esiny), 
~(—2E/a*)(n+-4ab/ R)}}, 


7 (—a/2E)(1 


From (6), 


-ecosy), 


nd t, is a constant. From (7), 


2( 2E \3 l—e|3 
Bi \?tan 1| J \2tan us 
al 1—e?| \1+e) = 
B 
= ji | F{B/x, sin—!(« cos q?)}— F{B/«, sin-1(a cos Q?)}], 
vhere y = 2!{— B/A—(B?/A®—4C/A)}} 


B = 24{— B/A+(B?/A? 


A 6ab R, 


B (n+-12ab/R), 
C n m? 6ab R, 


nd F(u,¢) denotes the incomplete elliptic integral of the first kind, viz. 








db 
: ; " dx 
F(u,¢) a ee Te 
J (l—p? sin’x)! 
: : , 0 
[The integral of (8), viz 
aq? 
; dx 
q3 a 
( . 9 > = * 9 
| {n—(m/sin x)?+ 6(ab/ R)sin®2}} sin®a 
Qi 


where Q5 


is a constant, can be expressed in terms of incomplete elliptic 
ntegrals of the third kind: these need not be given here as they are not 
required. 

Thus the motion is expressed in terms of six constants, viz. 

E, the total energy of unit mass, 

m, the angular momentum of unit mass about the earth’s axis, 

t,, the time of passage through apogee, 

Q?, the apogee distance 

Q2, the value of qg? for which g* vanishes, and 

n, which has in general no unique readily identified physical meaning 

but which enters most directly through (5). 

Naturally, some of these equations closely resemble the corresponding 
relations for Kepler motion. Thus e corresponds to the eccentricity of the 
Kepler ellipse, the apogee and perigee distances being (—a/2E)(1+-e) and 

-a/2E)(1 
radius if ( 


e) respectively. In particular, the orbit will have constant 


2E/a?)(n+-4ab/R) = 1, although such an orbit may not be 
confined to a plane. 


Kk 


5092. 40 
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3. The torsion 

Perhaps the most obvious effect caused by ellipticity of the figure of the 
earth is departure of the path of a freely falling body from a plane curve, 
although, of course, for certain obvious initial conditions the path will be 
plane. Turning aside from the special problem, it is a ready deduction from 
the Frenet reiations that the torsion of a free path in a field of acceleration 


ern 


is a 
J kr fi iam ii 


-_— (g™ Men | ae . 


where g”” is the (contravariant) metric tensor, 7” is unit tangent to the 





(9) 


path, e/* is the skew-symmetric (absolute) tensor whose non-vanishing 
} 


elements have magnitude g~!, g being the determinant |g,,,,,!, and /,..,, is 


the covariant derivative 
Cf, { n \ 
= 


| _— eq” a lie m 


For motion in a field of potential V, (9) is particularly convenient since 
f, = eV eq". If the field is symmetric about the axis of spherical polar 
coordinates, then 
ee ff T 
Kh fea —fefia)sin P+ (9)? fal fi ¢ sin g? +f, cos q?)}71 + 
Ufo2—fi 7) A sin q? —(fi,28in g?+-f, cos q?) fa}z?}7*. (10) 


9. 


An immediate consequence of (9) and (10) is that the torsion is zero if 


7? = 0 or if 71 = 7? = 0, i.e. if the instantaneous velocity is entirely in or 
normal to a plane containing the axis of symmetry. 
Applying (9) to the case of motion in the field of potential V,, and retaining 
only the first power of b, we find 
6bsin?q?{ x, t; | 


Ry r | (ay)?+ (ag)?) 


where x, denote the physical components of the tensor 7”. Hence |c! is 


(11) 


a maximum when 2, = 2,: in geographical terms, the projection of the path 
on the tangent plane to the earth’s surface (at the point having the same 
values of g? and q° as the particle) bisects the angle between adjacent pairs 
of cardinal directions on the compass. A similar application of (9) to the 
case of motion in the original field of potential V gives the result 
6b sin?q? (x, %3-+-2, 23 cot q?) 
(q)3 \ (x5)? i (a3)? : fe 





(12) 


4. Numerical estimates 
Finally, we may illustrate the influence of ellipticity of the earth’s figure 
by quoting some numerical values deduced from the relations given here. 
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Thus, suppose that a particle is set moving at right angles to the radius 
and in a meridian plane, with such speed that, in the Kepler field \j, a 
circular orbit would result. In the field V, it is found that 


é (2b/DR)|1—3 cos*0}, 


where D and © are, respectively, the values of q' and q? at the initial point. 
In particular, a circular orbit results if © = cos-!(-+1/V3), which corre- 
sponds to a geocentric latitude of 35° 20’ north or south. Although this 
orbit is circular it is not traversed at constant speed: it is found that the 
speed is (a/D)!{1+-(36/D R)sin?q?} and, if D and R are both of the order of 
the earth’s radius, the fluctuation in speed amounts to 0-16 per cent. On 
the other hand, the greatest value of e is found when © is 0 or z, and this 
value is about 0-0022 if R and D are of the order of the earth’s radius. 
Hence the total fluctuation in the distance of the particle from the centre 
of the earth amounts to about 28 km. This is, of course, much less than 
the expected variation in distance of an artificial satellite, due to the 
inaccuracy in controlling its initial velocity. 

For motion in the field of potential V,, and assuming the path to have the 
curvature of a terrestrial great circle, the ratio of the maximum magnitude 
f torsion to curvature is about 1/600. This may be otherwise expressed 
by giving the displacement that the particle undergoes in a direction normal 
to the osculating plane. If « is the curvature, the displacement parallel to 
the binormal is «os*/6, when the small displacement along the path is s. 
For s 1,000 km (roughly one-sixth of the earth’s radius) the displacement 
normal to the osculating plane is only 8 m. In calculating the effect of 
torsion in the path of a ballistic missile it is convenient to take the osculating 
plane at apogee as a reference; the displacements of the initial and final 
points of the path are equal in magnitude but oppositely directed normal 
to this plane. As an example we may consider a total path length of 
1.000 km, i.e. 2,000 km each side of apogee; the range on the surface of the 
earth will, of course, be somewhat less than 4,000 km. The curvature of the 
path at apogee is several times, perhaps four times, that of a terrestrial great 
ircle, and it is found that the maximum displacements of the initial and 
final points of the path, normal to the reference plane, have magnitude 
ibout 250 m; hence the relative displacement between those two points is 
ibout 500m. Thus, even in the most pronounced case, departure of the path 


from a plane curve is likely to be a matter of some refinement. 
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’ 


SUMMARY 
It is proved that Bessel functions of pure imaginary order (defined as a solution of 
equation (1) below) cannot have two zeros on a line parallel to the real axis. A 
consequence of this result is that Couette flow with static stability does not possess 
exponential unstable waves. It also follows from the proof that the Bessel function 
G,(z) defined by G,(z) 
the argument. 


hirie~*’™H\)(iz) can only have real zeros for finite values of 


1. Introduction 

BrEssEL functions of pure imaginary order play an important role in the 
study of some hydrodynamical models. An example is the following: the 
basic velocity is independent of the horizontal coordinates and a linear 
function of the vertical coordinate (Couette flow). The static stability is 
assumed constant, and the motion is limited by two horizontal planes. A 
discussion of the stability of this model has been given by Taylor (1) and 
recently by Eliassen, Héiland, and Riis (2). A result of their investigations 
is that in the case of positive static stability the problem of the existence 
of exponential unstable waves is reduced to the question of whether there 
exists a linear combination of two Bessel functions of pure imaginary order 
with two zeros on a line parallel to the real axis. The aim of this paper is to 
prove that such a combination does not exist, and the physical conclusion 
is that the model just described has no exponential unstable waves. 

If the uppermost plane (the vertical coordinate is positive upwards) is 
displaced so that the distance between the planes tends to infinity, the 
stability question depends on whether the Bessel function, of pure 
imaginary argument which vanishes for positive real arguments, has any 
complex zeros for finite values of the argument. It is shown below that this 
Bessel function has no such zeros, a result which has been published 
previously by Pélya (3) (only the result, not the proof, is published), in 
connexion with a discussion of the form of F(t) necessary to make 


@ 
| F(t)cos zt dt a function of z with only real zeros. 


0 
(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 4 (1957)] 
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ON THE ZEROS OF BESSEL FUNCTIONS 
2. Proof of the statement 
The equation for Bessel functions of pure imaginary order may be 


written . P 
d uldu, yim @ (1) 
dz? z dz \z? 


we 


where v isa real constant. A real value of z corresponds to Bessel functions 
of pure imaginary argument. As mentioned above we will prove that no 
solution of (1) can have two zeros on a line parallel to the real axis. To prove 
this we will make use of Green’s transform of a second-order differential 
equation. Consider a self-adjoint, linear differential equation of the 
second order, viz., d dw\ 


[K(2) )+ G(z)w = 0, (2) 


dz \ 


dz 


where K(z) and G(z) are assumed to be analytical functions in a domain D 
throughout which K(z) does not vanish. If 


dw 


w w, = K(z) ie” 
dz 


1 ? 2 


the Green’s transform of equation (1) is (see (4)): 
[ow 7 — | |w,|*L(z) dz | w, |*G(z) dz = 0, (3) 


assuming that every point of the path of integration lies in D. Here a bar 
denotes the complex conjugate of the quantity, and L(z) is defined by 


L(z) = 1/{K(z)}. 


Separating the real and imaginary parts, we have 


Ze2 <2 
re[@, w.|*—re | |w,|L(z) dz+-re | |w,|?G(z) dz = 0, (4) 


ahs) 


n 
— 
~ 
x 
i! 
| 
~ 
—) 
— 
| 
~_ 


im|[@, wa)?—im | |w,|2L(z) dZ+im | |w,|2G( 


Ne 


We first prove that, if a solution of (1) with two zeros on a line parallel 
to the real axis exists, then the two zeros must lie on opposite sides of the 
imaginary axis (it will be assumed throughout the paper that the two zeros 
can be combined with a path which does not enclose the origin). To prove 
this statement we assume for the moment that the two zeros are on the 
same side of the imaginary axis. By introducing 


u zty 
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equation (1) takes the self adjoint form 


d dv 
a2 - 2__ 42\y; - 0 a 
dz (: i er ” 


> 


where py? = v4 


pp 


Comparing (7) with (2) we notice that in our case 


| l x2—y?— Jay 


i ha 9 > 9\° > 
mi2Z) 2 (a?+-y?)? 
G(z) p2—z2? p?—a?+ y?—2aryt, 
since Z2= x+y. 
Let z, and z, represent the two zeros. In (5) the first term vanishes. As the 
path of integration we choose the straight line connecting z, and z,. It is 
then seen that the two integrals in (5) have opposite signs, and thus (5) 
cannot be satisfied. 
To complete the proof we return to (1) and now introduce 


u 2 *Ww. (9) 
Equation (1) takes the form 


4 (fe Jw =— (0). (10) 


This equation is self adjoint with 


2 2 42 __ Doras 
and Q(z) =F .—| pe" y manatee —1. (11) 
Z (x*-+-¥*)* 


We assume again that a solution with two zeros on a line parallel to the real 





's exists. From the above result these two zeros must lie on opposite 
sides of the imaginary axis. Choosing z, and z, to represent the zeros, the 
first term in (4) vanishes. There are two possibilities for the path of integra- 
tion: (i) if the imaginary parts of the zeros are positive, we integrate along a 
line parallel to the imaginary axis in a positive direction to a point A, then 
along the line parallel to the real axis to a point B, and finally along a line 
parallel to the imaginary axis to the other zero point; (ii) if the imaginary 
parts of the zeros are negative, the path of the integration will be the 
corresponding path in the half-plane of y < 0. Let z, represent the zero 
with negative real part. It is then readily seen that the second term in (4), 


namely, —re 





w, *L(z) dz, will be either zero or negative on the path of 
zy 


integration. The third term, re | |w,\?G(z) dz, is negative on the lines 


parallel to the imaginary axis. By choosing the distance between the line 
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ON THE ZEROS OF 





BESSEL 





through the zeros and the line AB large enough, it can be seen that the 
contribution from the path AB is also negative. Equation (4) cannot be 
satisfied, and we conclude that two such zeros cannot exist. 

The arguments given above are also valid if one of the zeros tends to 
infinity. We therefore see that the function which is often denoted by 
@(z) and defined by 

G,(Z) bare tua FD ({z), (12) 
where H‘» is a Bessel function of third order, has no complex (or pure 
imaginary) zeros for finite values of the argument since this function tends 


to zero when z — ooand /arg z }7. Asimilar result is true for the solution 
of (1) which tends towards zero when z > 0 and $a < argz < 37. 


It should perhaps be mentioned that these functions have an infinite 
number of real zeros. 

By applying a similar method to that used above it follows also that no 
solution of (1) can have two zeros which lie on a line parallel to the imaginary 
xis and on the same side of the real axis. To obtain the results above we 
have made use of the fact that 

Be ie 

isa positive quantity. The results derived are also valid for pure imaginary 


values of v, provided that . 1 
yp” . 4° 
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DETERMINATION OF THE OPTIMUM RESPONSE OF 
LINEAR SYSTEMS 
(ZERO DISPLACEMENT ERROR SYSTEMS)+ 
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SUMMARY 


In (1) two response functions Z and L, were defined by the equations 


x oO 
r (de\? 
L e?dr and L, = | ( : dr, , 
P J \dr : 
0 
where ¢ is the error at time 7. These relations are applied to a linear system having 


one degree of freedom. By considering the response to a step-function disturbance 
it is found that systems making LZ a minimum have a lightly damped oscillatory 
response. The smaller L, is, the ‘smoother’ is the response. 

Values are obtained for the coefficients of the characteristic equation of any order 
making J a minimum. An approximate method is given for correcting these coefii- 
cients to enable the response to be improved to give equal damping in the least damped 
modes of oscillation. 


1. Introduction 
In (1) two response functions L and L, were defined by the equations 


x 


L= [ e? dr (1) 
0 
and L, = | (5) ar, (2) 
pe at 


0 

where e is the error at time 7. Formulae were given for L and L, 
in terms of (i) the coefficients of the characteristic equation, and (ii) the 
frequency response spectrum. We shall now show how L and L, can be used 
to obtain the optimum response. As stated in (1) a system based on 
minimizing L often has an undesirable overshoot. We shall consider the 
performance of systems for which L is near its minimum value while J, is 
considerably lower than its value when L is a minimum. We shall obtain 
the relation between the transient response, the frequency response, and 
the roots of the characteristic equation for systems having no static error. 

We consider a system for which the equation of motion is 

n n—-ly N—2y 

ay, a a be +vag + ve Fy = + Agt = f (7), (3) 
+ The material in this paper is taken from a Ph.D. thesis for the University of Glasgow. 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 4 (1957)] 
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where dp, @, @,..., a, are constants and a, + 0, and z is the response of 


n 


the system with input f. 


2. Response of zero displacement error systems to step function 
disturbance (constant displacement input) 
In this case 
f=0 for7 <9; f=f, forr> 0, (4) 
where f, is a constant. 
As in (2) we shall find it convenient to normalize (3). 
Let Ay = wHa,, (5) 
and define a new time scale by the relation 
u WoT. (6) 


From (3) to (6) the normalized equation of motion is 


d”x d"-ly d”—*x dx ho = 
T Un-1 i 1 4n-2 gee TH Te (r > 0), (7) 

du’ du" du"-* du Ay 
where q, = a,/(a, wa") (r= Oton). (8) 


Now the magnitude of the response x is proportional to f,, but the nature 
of the response (per cent overshoot, damping time) is independent of f,. 


We shs ake , 
e shall take fy ay. (9) 


Thus we shall determine the response of the system given by (7) to a unit 
step function disturbance. We shall only consider stable systems; then the 
static error ultimately vanishes since x > 1. Such systems are sometimes 
called ‘zero displacement error systems’ (see (2)). As shown in (1) the 
response is identical in form with that in a free motion with x, = —1 and 
Dz, D*z,,..., D®—'z 
7 = 0, and D is the differential operator. 


» zero, where the suffix 0 denotes the values at time 


3. Second-order zero displacement error system 


From (1), with 2, 1, Dx = 9, 


2La,a, a, a, az+-dy Qs (10) 
ay ay 


or, in terms of the normalized coefficient g, given by (8), 


214 
2 Lw, . -. (11) 
1 
; 2L ] 
Similarly, - Ass =—. (12) 
Wy 11 


Now for stability with a, > 0, both a, and a, must be positive, i.e. gq, > 0. 
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We see that for a system with a given undamped natural frequency w,. 

g 1 J % 

L and L, are both functions of g, (which for a mechanical system is pro- 
1 11 : I 


portional to the damping). The minimum value of L occurs for q, = 1, ie, 


Gt = gt. 
Then = L/w», (13) 
and for this value of qg, L, = 4up. (14) 


This corresponds to an oscillatory motion with overswing of 16 per cent 
(see Fig. 1). We shall see that this tendency for the oscillatory motion to 


Sr 20; 1200 





° 
: oa cy) 
* 1-0F < o = 
g S 4 
c 
S : : 
2 O57 - 0 |-200 2 
. + © “40705 10 20 50° 
Non dimensional time a7 Non dimensional frequency 


Fic. 1. Second-order zero displacement error system 


be rather lightly damped is a characteristic of systems based on L This 
overswing is lessened by selecting a higher value of q, (or a,) as is seen from 
the following table. 


min* 


TABLE | 


Variation of response functions and overshoot with the damping 
coefficient for a second-order zero displacement error system 





(dp Ag) Lewy L, Wo °%, overshoot 
6 I'l o'82 28 
5 5 3 
o'5 1°03 0°63 26 
I*o 1°00 0°50 16 
I*2 1°02 0°42 9 
I°4 1°06 o* 30 5 
I°O I°it o°3I 2 
1°8 1°18 0°25 fe) 
C 1°25 "25 ° 





The minimum of L is seen to be fairly flat; the value of L is not more 
than 10 per cent greater than the minimum for values of g, from 0-7 to 1-5. 
We see that there is good correlation between the per cent overshoot and 
the value of L,, as would be expected from the definition of L,. The smaller 
L, is, the ‘smoother’ is the response, which becomes less oscillatory in 
nature. Precisely what value of q, is chosen for the most satisfactory system 
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sa matter for individual choice, depending upon the acceptable per cent 
wershoot. From the above table we see that to keep both L and L, as 
small as possible, g, should be greater than 1. Now if q, is greater than 2, 
the roots of the characteristic equation will both be real and negative, one 
root decreasing, the other increasing, as q, increases; the motion will then 
be composed of two monotonic functions (subsidences) and, in the termino- 
ogy of servomechanism analysis, would be considered to be overdamped. 
We are therefore led to the criterion 

lO<q, < 2°0 (15) 
for satisfactory performance for a second-order zero displacement error 
system. This is in good agreement with current practice with servo- 
mechanisms. 

The choice of w, will depend on the desired speed of response of the 
system; the higher the value of wy. the sooner will the system reach its 
steady state. Changing w, merely alters the time scale of the damping but 
loes not affect the per cent overshoot or the general form of the response. 

The other diagram in Fig. 1 shows the gain (in decibels) plotted on a 
garithmic scale against the non-dimensional frequency w/w, for the two 
seS 4, | and q, = 2. If M is the dynamic magnification (1), the 
gain 20 log,) VM decibels. 

For the range of q, given by (15), M,,,, does not exceed 1-15, which is 


max 


well within current practice with servomechanisms. 


4. Higher-order zero displacement error systems 


The above analysis can easily be extended to higher-order systems. From 


1), with Xo - 
Dx, = D*x, = ... = D"-1z, = 0, 
2fa,| a, ds Gg Gy .. . | @, @y @& G& . . .« | (16) 
A A, a, a6 ' ‘ R ar) a as as : : 
0 @ @ @ .. . O Gy. Ge 
0 Ay Ag UN . . ' 0 0 a, as, 


or, in terms of the normalized coefficients q,, qo,.--, ¢n-, given above, 


2Lwo| 9, 93 Ws 92: | |e @& & & +: « . | (17) 


I q< Ga 6 = I 1 qs W5 . . ° 
LU 11 Ws 45 co 8 oe | 0 I G2 GW <:-+ «+ - 
® bt @ @ « « «wf . &©§ & @& « « «3 


where the determinants on the left-hand sides of (16) and (17) are of order 
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n—1 and those on the right-hand sides are of order n. Also q,, = 
Similarly, 


21, % %% % G% + - -|=| GW % Ww WW - + + | (18 
wy| | Ms Ge GW + «+ * | % Ws Ws da + - 

9 & WW IG + - - 1 Yo Us 6 . 

O 1 G@ Ww + + - | 9 % % %W +: - 


| . . . . . . . . | . . . . . . 


the final determinant on the right-hand side of (18) being of order n—3, 
Equation (18) can be written in a simple form in terms of the test functions 
of the characteristic equation. We find 
7 
I = QQ T-2 (19 
“2 27 ? . 
oen—l 
where 7, is the rth test function of the characteristic equation corresponding 
to (3) (see (3) and (4)). For stability when a, > 0, all the test functions 
must be positive. 
The minimum values of L and the corresponding values of qo, 41, Yos-++ 4, 

are given in the following table. 


TABLE 2 


Values of L,,;,, and the normalized coefficients q,, for zero displacement 
error systems of order n 


n 2L nino do nN 
I I I I qe 
2 2 I I I q3 
3 3 I 2 I I #6. cy 
4 4 I 2 3 I I qs 
5 5 I 4 3 4 I I I6 
6 6 I 3 6 $ 5 I I q7 
7 7 I 4 6 10 5 6 I I 18 
8 8 I 4 10 10 15 6 7 I 1 
From the above table we see that 
n 
9 
} (20 
-Wo 


and that the.coefficient q,, for an nth order system is equal to 


p 
binomial coefficient, namely 


Y ! ! n 
vm = p!/{m!(p—m)!} 








m+n . cars 
asanues if m+n is even 
5 (21 
and Pp = 
m+t+n—l . , 
—;—_ if m+n is odd 


Thus the above table can easily be extended to any order n. 
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We find also that when the test functions are written in terms of the q’s 


min? 
7, = apport §=(¢ = I...., 8). (22) 
From (19) and (22) we find that when L is a minimum 
] ai XM Wo, (23) 
“1 9,,,n lq > a 
-Wo n - 


min 18 proportional to the order n of the system, but the corresponding 
value of L, is independent of the order of the system. 

The motion for this system is composed of a number of damped oscilla- 
tions (and possibly a monotonic function if m is odd), as shown in the 
following table. 

TABLE 3 
Roots of the characteristic equation for a zero 


displacement error I system of order n 


“min 
Root 
7 IS+1°312 
5 505 o°105+-1°571 
} 235+-0°Sd1 0°06 -+- 1°701 
5 62 O°155-1°502 
1°78 
0°22+-0°06051 0°09 + 1°35 
5 1°52 
27 +0°2531, —O°I5+0°9I1, 
63+ 1°501, —0°013+1°501 


For a system of given order the modes of low frequencies are the more 
highly damped. As n increases, the damping in all the modes decreases 
and the corresponding frequencies increase. The mode with the least 
damping is an oscillatory mode of relatively short period. As would be 
expected, this gives a system with rather a wavy response, due to the 
terms of high frequency, with a certain amount of hunting (usually of 
small amplitude compared with the input motion). This is shown in Fig. 2 
for systems of fourth and fifth order, the first overswing for a unit step 
disturbance being 9 per cent for the fourth-order system and 10 per cent 
for the fifth-order one. The amplitudes of the corresponding motions at, 
SAY, WoT 10 are 6 and 11 per cent respectively. 

As with the second-order system it is possible to find a more satisfactory 
_ by choosing the q’s so that L, is decreased 
This could be done by 


response than that for L,,,;, 
considerably while L is only increased slightly. 
const. 


determining the locus at which the surfaces L = const. and L, = 


touch one another. However, we shall be mainly interested in points close 


to L 


min 


and the following approximate procedure can be adopted. 
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For values of the q’s which make J a minimum we find that 


: 1 —9 (r = l,..., n—3) 
eq, 
. (24) 
OL | 
: 1 — —uy (rj n—2,n—1) 
eq, 


Thus for points in the neighbourhood of L,,;,, £2, can be most effectively 
reduced by increasing q,,_, and q, _, by some quantity, @ say, i.e. by choosing 
values of the q’s which are on the normal to the surface L, = 0-5w, at the 


point corresponding to L,,in- 






ie 
| id 
§ « 8 
2 2 
c | c 
a | 4 
BC sit | 3s 





0 5 0 5 0% 5 
Non dimensional time a7 








Non dimensional time wo7 


Fic. 2. Fourth and fifth order.zero displacement error system 


As @ increases, the roots of the characteristic equation vary, the damping 
of the lightly damped oscillation being increased while that of the other 
modes is either unchanged or slightly decreased. We note that since q,,_, 
increases as @ increases the ‘total’ damping of the system will increase. For 
a particular value of # the characteristic equation has at least three (usually 
four) roots with equal negative real parts. For such a system L is denoted 
by Ly. The values of the normalized coefficients and the roots of the 
characteristic equation are given in Tables 4 and 5, obtained by the approxi- 
mate procedure outlined above. 


TABLE 4 


Values of the normalized coefficients q,, for zero displacement error systems 


of order n with equal damping (Ly system) 


n Yo bt 12 q3 

3 I 2°5 I°5 I q4 

4 I 2 3°35 1°38 I qs 

5 I 3 3 $°17 1°17 I I6 

) I 3 6 4 5°07 1°07 I q7 

7 I 4 6 | fe) 5 6°03 1°03 I Is 
5 I j 10 10 5 6 7°O14 I°O14 I 


We see that, as the order of the system increases, the value of 6 for equal 
damping decreases. Comparing Tables 3 and 5 we see that in passing from 
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Li, to La the damping of the least damped oscillation is increased about 
three times, the frequencies of the various modes being practically un- 
changed. Thus as shown in Fig. 2 the response of the system to a unit step 
disturbance is of amuch smoother nature. The peak overshoot has increased 


slightly to 13 per cent for both systems. 
TABLE 5 


Roots of the characteristic equation for a zero 


displacement error La system of order n 


, Root 
5 I*32 
} 53 345 = 1°54 
I ) I9+1°7 
$7 1O5-+-1°172 
I°77? 
22 7Ol 0°005 + 1°3601 
1°S3] 
” 282 C 


Comparing the response curves of Fig. 2 with those based on the I.T.A.E. 
criterion (2) and on the Butterworth filters, (5) and (6), we see that the 
vershoot is greater, and the damping less, than with the I.T.A.E. systems. 
The overshoot is about the same as with the Butterworth filters, the 
damping is slightly less. The time for the error first to become zero 
momentarily) is smaller than for either the I.T.A.E. systems or the 
Butterworth filters (for a given value of wo). 

The advantage of the above determination over that of the I.T.A.E. 
method is that the coefficients given in Table 4 can be readily extended to 
systems of any order, being based on an analytical formula. The response 
curves of the systems L, and L,,;,, form a family of curves as the order of 
the system is increased. The justification of the use of the L criteria must 
depend on the final form of the response curves and upon the comparative 
ease with which the criteria can be calculated, especially when the approxi- 
mate method is used. It is not surprising that the general shape of some 
of the response curves for L,, is not unlike that with the Butterworth filters. 
Other investigators have arrived at similar results (see the discussion of (2)). 
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CORRIGENDUM 
‘A SOURCE IN A ROTATING FLUID’ by 8S. N. BARUA 


In the paper mentioned above (Quart. J. Mech. App. Math. 8 (1955) 27) the sentence 
in which equation (11) appears should read as follows: Since J,(a, c/a) tends to a, ¢/2a 
(instead of «a, c/a) as c tends to zero, the expression for d may be written as 
W % l 
d . ~« > ——. -e-*n@/A J (y 7 a). 
27a? Lot Xn J? (Xp) hee 
Consequently the equation (12) of the paper becomes 


| 2 fa 2 Lqd Xn 0/a ‘Sis Xn 0/A\ 2 
(“= “) (ete) stog(2t*) 0-3918| 2 Poe So | 

a \ a . a —_—sSd* (x, Lat J (x,)/ J’ 
and Table 1 becomes 


° 
4 0°2845 o's 
I 0*5106 X 107} 0°28 





These results are illustrated in Fig. 1. 
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Fie. 1. A section of the discontinuity-surface. (The plotted points on the 
discontinuity-surface are the results of Table 1) 
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CORRIGENDUM 


This figure was omitted from * Experiments on the Slow Swirling Flows 
of a Viscous Liquid Through a Tube’ by A. M. Binnie, Part 3 of 
Volume X 
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and Physicist. Applied Analysis is a philosophical and strictly 
theoretical approach to mathematical methods used in the solution 
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